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1.  Introduction 

In  this  paper,  we  develop  and  analyze  robust  multilevel  preconditioners  for  discontinuous 
Galerkin  (DG)  discretizations  of  the  second  order  elliptic  equation  with  strongly  discontinuous 
coefficients: 

/  -V  ■  (kVu)  =  f  in  n,  n 

\  u  =  0  on  dQ,  ' 

where  hi  C  IRc/  is  a  bounded  polygon  (for  d  =  2)  or  polyhedron  (for  d  =  3).  The  scalar 
function  n  =  k{x)  denotes  the  diffusion  coefficient  which  is  assumed  to  be  piecewise  constant 
with  respect  to  an  initial  non-overlapping  subdomain  partition  of  the  domain  Q,  denoted  Ts  = 

_  _  O  O 

{f with  Um=i^m  =  G  and  n  Dn=  0.  Although  the  (polygonal  or  polyhedral)  regions 
VLm  ,m  =  1 ...  M  might  have  complicated  geometry,  we  will  always  assume  that  there  is  an  initial 
triangulation  7o  such  that  kt  =  k(x)\t  is  a  constant  for  all  T  €  %■  Problem  (1.1)  belongs  to 
the  class  of  interface  or  transmission  problems,  which  are  relevant  to  many  applications  such  as 
groundwater  flow  [54],  electromagnetics  [53],  semiconductor  device  modeling  [31,  59],  and  fuel 
cells  [74,  76].  The  coefficients  in  these  applications  might  have  large  discontinuities  across  the 
interfaces  between  different  regions  with  different  material  properties. 

Finite  element  discretizations  of  (1.1)  lead  to  linear  systems  with  badly  conditioned  stiffness 
matrices;  not  only  is  the  ill-conditioning  with  respect  to  the  mesh  size,  but  the  condition  num¬ 
bers  depend  linearly  on  the  largest  jump  in  the  coefficients.  Much  research  has  been  devoted 
to  developing  efficient  and  robust  preconditioners  for  use  with  iterative  methods  for  conforming 
finite  element  discretizations  of  (1.1).  The  resulting  preconditioners  can  be  grouped  into  two 
basic  classes  of  methods:  non-overlapping  and  overlapping  methods.  Both  classes  of  precondi¬ 
tioners  have  their  advantages  and  disadvantages,  and  as  a  result  it  is  desirable  to  have  access 
to  robust  and  efficient  methods  from  both  classes.  Domain  decomposition  Balancing  Neumann- 
Neumann  [58],  FETI-DP  [55]  and  Bramble-Pasciak-Schatz  preconditioners  [12]  belong  to  the 
first  class  of  methods  (the  non-overlapping  methods).  They  have  been  shown  to  be  robust  with 
respect  to  coefficient  variations  and  mesh  size  (up  to  a  logarithmic  factor),  in  theory  as  well  as 
practice,  but  only  if  special  coarse  spaces  and  coarse  solvers  (such  as  those  based  on  discrete 
harmonic  extensions  [42,  58,  43,  65])  are  constructed  (see  also  a  survey  paper  by  Xu  and  Zou  [79] 
for  many  techniques  and  references). 

The  class  of  overlapping  methods  encompasses,  among  others,  overlapping  Schwarz  meth¬ 
ods,  geometric  multigrid,  and  more  general  Multivel  methods,  as  such  the  Bramble-Pasciak-Xu 
(BPX)  preconditioner  (see  e.g.  [13,  14,  77]).  In  practice,  it  has  always  been  observed  that  all 
these  methods,  when  used  as  preconditioners  in  conjugate  gradient  iteration,  result  in  efficient 
algorithms  that  behave  robustly  with  respect  to  the  jump  in  the  coefficients,  independently  of 
the  problem  dimension.  However,  the  first  analyses  show  robustness  only  in  certain  particular 
cases:  quasi-monotonicity  assumption  on  the  coefficients  [41];  space  dimension  d  <  2  [10];  or 
when  each  subdomain  Q,m  touches  the  Dirichlet  boundary  [75,  61].  In  general,  these  theories 
predict  a  deterioration  in  the  rate  of  convergence  of  multigrid  and  overlapping  domain  decom¬ 
position  methods,  with  respect  to  both  the  coefficients  and  the  mesh  size.  An  improvement  in 
the  rate  of  convergence  can  be  achieved  by  resorting  to  basis  stabilization  techniques  [70],  or  by 
employing  energy  minimizing  coarse  spaces  [73]. 

However,  we  note  that  there  is  a  discrepancy  between  these  theoretical  results  and  the  con¬ 
vergence  rate  actually  observed  in  practice.  This  is  due  to  the  fact  that  the  analysis  approach 
used  in  earlier  works  was  based  mainly  on  the  standard  theory  of  CG  (see  for  example  [51,  The¬ 
orem  9.4.12],  or  [63,  Theorem  6.29])  which  provides  only  bounds  in  the  worst  case  scenario,  and 
does  not  exploit  the  spectral  structure  of  particular  preconditioned  system.  The  preconditioned 
system  may  actually  have  only  a  few  small  eigenvalues  (depending  on  the  coefficient  distribution 
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and  mesh  size),  while  the  other  eigenvalues  are  bounded  nearly  uniformly.  This  was  first  ob¬ 
served  by  [48,  72]  for  the  simple  Jacobi  preconditioner.  By  using  a  more  sophisticated  approach 
to  CG  theory  (see  [6,  Section  13.2],  [7])  and  considering  the  distribution  of  the  spectrum  of  the 
preconditioned  system,  it  is  possible  to  show  that  the  small  eigenvalues  do  not  influence  the 
(observed)  asymptotic  convergence  rate.  This  approach  has  been  pursued  in  [78,  80],  where  it 
is  show  that  both  standard  multilevel  and  overlapping  domain  decomposition  methods  lead  to 
nearly  optimal  preconditioners  for  CG  algorithms.  See  also  [26]  for  further  extensions  in  the 
case  of  highly  graded  adaptive  meshes. 

While  preconditioners  for  conforming  discretizations  have  been  considered  by  a  number  of 
researchers,  less  work  has  focused  on  the  construction  of  preconditioners  for  non-conforming 
approximations  of  (1.1).  In  this  article,  we  consider  preconditioners  when  the  family  of  Interior 
Penalty  (IP)  Discontinuous  Galerkin  (DG)  methods  are  used  to  approximate  (1.1).  DG  methods 
are  designed  to  provide  robust  discretizations  of  partial  differential  equations  of  nearly  any  type, 
and  can  even  handle  equations  whose  type  changes  within  the  computational  domain.  They  are 
naturally  suited  for  multi-physics  applications,  and  for  problems  with  highly  varying  material 
properties,  such  as  problem  (1.1).  The  development  of  efficient  solvers  for  DG  discretizations  has 
been  pursued  only  in  the  last  ten  years;  however,  due  to  the  growing  importance  of  DG  methods, 
there  is  now  considerable  interest  in  this  area.  Examples  of  overlapping  preconditioners  for  non- 
conforming  discretizations  are  found  in  [65,  64],  where  the  analysis  depends  on  the  assumption 
that  the  coefficient  k  is  quasi-monotone.  Since  the  algorithms  we  develop  here  contain  the 
solution  of  systems  corresponding  to  non-conforming  discretizations  of  the  model  problem  (1.1) 
as  a  part  of  the  preconditioning  step,  the  analysis  presented  here  gives  bounds  on  the  convergence 
rate  for  two  level  and  multilevel  methods  for  the  lowest  order  non-conforming  Crouseix-Raviart 
finite  element  discretizations. 

While  classical  approaches  have  been  successfully  extended  for  second  order  elliptic  prob¬ 
lems,  discontinuous  nature  of  the  underlying  finite  element  spaces  has  motivated  the  creation 
of  new  techniques  for  the  design  of  solvers.  Additive  Schwarz  methods  (of  overlapping  and  non¬ 
overlapping  type)  are  considered  and  analyzed  in  [44,  38,  2,  3,  4,  11],  Multigrid  methods  are 
studied  in  [47,  20,  19,  18,  62,  30].  Two  level  methods  are  presented  in  [34,  22,  23].  More  general 
multi-level  methods  based  on  algebraic  techniques  are  considered  in  [57,  56].  However,  except 
from  several  numerical  experiments  reported  in  [34,  33,  3,  35,  56],  all  these  works  deal  with  the 
case  of  a  smoothly  or  slowly  varying  diffusivity  coefficient.  For  problem  (1.1),  only  in  [38,  39,  40] 
the  authors  introduce  and  analyze  non-overlapping  BBDC  and  FETI-DP  domain  decomposition 
preconditioners  for  a  Symmetric  Interior  Penalty  discretization  of  (1.1).  The  DG  discretization 
in  only  used  on  the  skeleton  of  the  subdomain  partition,  while  a  standard  conforming  approx¬ 
imation  is  used  in  the  interior  of  the  subdomains.  Robustness  and  quasi-optinrality  is  shown 
for  the  Additive  and  Hybrid  BBDC  [39]  and  FETI-DP  [40]  preconditioners,  even  for  the  case 
of  non- matching  grids.  The  analysis  of  the  first  preconditioner  requires  an  interface  condition 
relating  the  magnitude  of  the  coefficient  and  the  mesh-size. 

The  goal  of  this  article  is  to  design,  and  provide  a  rigorous  analysis  of,  a  simple  multilevel 
solver  (belonging  to  the  overlapping  category),  for  the  family  of  Interior  Penalty  methods  for 
approximating  (1.1).  Together  with  the  family  of  IPDG  methods  (including  symmetric  and 
non-symmetric  methods),  we  consider  the  corresponding  family  of  weakly  penalized  or  IPDG-0 
methods  (called  Type-0  in  [9] ) .  Our  approach  follows  the  ideas  in  [9] ,  and  it  is  based  on  a  splitting 
of  the  DG  space  into  two  components  that  are  orthogonal  in  an  energy  inner  product  (more 
precisely,  in  the  energy  inner  product  induced  by  the  IPDG-0  methods).  Roughly  speaking, 
the  construction  amounts  to  identifying  a  “low  frequency”  space  (Crouzeix-Raviart  elements) 
and  then  defining  a  second  complementary  space.  Such  a  decomposition  turns  out  to  be  critical 
to  the  design  of  robust  preconditioners  in  the  case  of  PDEs  with  rough  coefficients,  such  as 
problem  (1.1).  However,  a  notable  difference  takes  place  in  the  decomposition  of  the  DG  space 
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introduced  for  the  Laplace  equation  [24,  9].  For  problem  (1.1),  a  stable  decomposition  can  be 
obtained  by  introducing  subspaces  that  depend  on  the  coefficient  n,  and  this  is  certainly  related 
to  the  splittings  used  in  algebraic  multigrid  (AMG  [16]),  where  one  seeks  space  decompositions 
depending  on  the  operator  at  hand.  In  addition,  both  components  of  the  splitting  have  a  locally 
supported  basis.  As  we  will  show,  the  analytical  results  and  the  preconditioning  techniques 
obtained  through  the  decomposition  introduced  in  [9]  for  the  Laplace  equation  can  be  extended 
to  the  case  of  problem  (1.1)  with  jumps  in  the  coefficients.  With  the  orthogonal  splitting  of 
the  DG  space  at  hand,  the  solution  of  problem  (1.1)  reduces  to  solving  two  sub-problems;  a 
non-conforming  approximation  to  (1.1),  and  a  problem  in  the  complementary  space  containing 
high  oscillatory  error  components.  We  show  that  the  latter  problem  is  easy  to  solve,  since  it 
is  spectrally  equivalent  to  its  diagonal  form,  and  therefore  CG  with  a  diagonal  preconditioner 
is  a  uniform  and  robust  solver.  For  the  former  problem,  that  is  the  approximation  in  the  low 
frequency  or  Crouziex-Raviart  space,  we  develop  and  analyze  a  two-level  (overlapping)  method 
and  a  BPX  preconditioner.  In  fact,  the  theory  for  the  BPX  type  multilevel  preconditioner 
follows  easily  from  the  analysis  of  the  two-level  method. 

We  follow  the  approach  taken  in  [78,  80]  involving  estimating  the  asymptotic  rate  of  con¬ 
vergence  of  the  preconditioned  system.  Nevertheless,  dealing  simultaneously  with  the  jump  in 
the  coefficient  k  and  the  non-nested  character  of  the  CR  spaces  presents  extra  difficulties  in 
the  analysis  which  preclude  from  a  simple  extension  of  the  works  [78,  80].  We  are  able  to  es¬ 
tablish  nearly  optimal  convergence  and  robustness  (with  respect  to  both  the  mesh  size  and  the 
coefficient  k)  for  the  two-level  method  and  for  the  BPX  preconditioner  (up  to  a  logarithmic 
term  depending  on  the  mesh  size).  The  resulting  algorithms  involve  the  use  of  a  solver  in  the 
CR  space  that  is  reduced  to  a  smoothing  step  followed  by  conforming  solver.  Therefore,  in 
particular  one  can  argue  that  any  of  the  robust  and  efficient  solvers  designed  for  conforming 
approximation  of  problem  (1.1)  could  be  used  as  a  preconditioner  here.  The  solvers  are  shown 
to  be  robust  and  uniformly  convergent  for  symmetric  IPDG-1  methods.  For  the  non-symmetric 
IIPG  and  NIPG  we  propose  two  preconditioners;  the  symmetric  part  and  a  block- Jacobi- two 
level  method.  For  the  former,  we  indicate  how  the  theory  for  second  order  problems  in  [9]  can  be 
successfully  extended  to  cover  problem  (1.1),  once  the  provisions  given  in  this  paper  are  taken 
into  account.  Unfortunately,  for  the  latter,  in  spite  of  its  simplicity,  we  are  not  able  to  provide 
a  complete  theory  at  this  time,  and  it  will  be  subject  of  future  research. 

Outline  of  the  paper.  The  rest  of  the  paper  is  organized  as  follows.  We  introduce  the  IPDG-1 
and  IPDG-0  methods  for  approximating  (1.1)  in  §2  and  revise  some  of  their  properties.  The 
space  decomposition  of  DG  finite  element  space  is  introduced  in  §3.  Consequences  of  the  space 
splitting  and  solvers  for  the  IPDG-0  methods  are  described  in  §4.  The  two-level  and  multi-level 
methods  for  the  Crouziex-Raviart  approximation  are  constructed  and  analyzed  in  Section5.2 
(see  §5.3).  Section  6  is  devoted  to  the  design  and  analysis  of  the  solvers  for  both  the  IPDG-0 
and  IPDG-1  family.  Numerical  experiments  are  included  in  §7,  to  verify  the  theory  and  assess 
the  performance  and  robustness  of  the  proposed  preconditioners.  The  paper  is  completed  with 
an  Appendix  where  we  have  collected  proofs  of  several  technical  results  required  in  our  analyses. 

2.  Discontinuous  Galerkin  Methods 

In  this  section,  we  introduce  the  basic  notation  and  describe  the  DG  methods  we  consider  for 
approximating  problem  (1.1).  To  begin,  given  a  triangulation  7/  of  the  domain  f 2,  we  denote 
£h  :=  £f  U  the  set  of  all  edges  (2D)  /faces  (3D)  of  7/,  where  £%  is  the  set  of  all  interior 
edges/faces,  and  £ ®  is  the  set  of  all  boundary  edges/faces.  Throughout  the  paper  we  shall  use 
the  standard  notation  for  Sobolev  spaces  and  their  norms.  We  denote  H2(Th)  as  the  set  of 
element-wise  H2  functions,  and  denote  L2(£h)  as  the  set  of  L 2  functions  defined  on  £)x.  We  will 
also  use  the  notation  x\  <  y\ ,  and  X2>i/2,  whenever  there  exist  constants  C\ ,  C2  independent 


MG  FOR  DG  WITH  JUMP 


5 


of  the  mesh  size  h  and  the  coefficient  k  or  other  parameters  that  aq,  xq,  y i  and  y2  may  depend 
on,  and  such  that  x\  <  C\y\  and  xq  >  C^yi- 


Trace  Operators.  Following  [5] ,  we  recall  the  dehnition  of  the  average  and  jump  trace  opera¬ 
tors  for  scalar  and  vector-valued  functions.  Let  T+  and  T~  be  two  neighboring  elements,  and 
n+,  be  their  outward  normal  unit  vectors,  respectively  (n^  =  nr±).  Let  ^  and  t±  be  the 
restriction  of  £  and  r  to  T±.  We  set: 

ia  =  \((+  +  C),  [Cl  =  C+n+  +  rn”  onee£°h,  (2.1) 

M  =  \{T+  +  r“),  [t]  =  r+  •  n+  +  T~  •  n"  on  e  G  8%, 

We  also  define  the  weighted  average,  -{[•]}■  5,  for  any  6  =  {<5e}ee££>  with  5e  G  [0, 1]  Ve,  : 

{{Cl<5  =  5eC+  +  (1  -  $e)C  ,  =  5er+  +  (1  -  5e)r~  ,  on  e  €  £%  .  (2.2) 

For  e  G  £®,  we  set 

[Cl  =  Cn,  f  rj-  =  =  r  on  e  G  £%.  (2.3) 

We  will  also  use  the  notation 

(u,  w)Th  =  E  /  uwdx  V  u,  w  G  L2(fi),  (u,w)£h  =  ^2  uw  V«,ro,6l/2(4). 

TeTh  Jt  ee£h  Je 

The  DG  approximation  to  the  model  problem  (1.1)  can  be  written  as 

Find  Uh  G  VDG  such  that  Adg(u,vj)  =  (f,w)j-h  ,  Vic  G  VDG  ,  (2.4) 

where  ADG(-,  •)  is  the  bilinear  form  defining  the  method  and  can  be  of  two  different  types.  For 
the  Type-1  family  of  weighted  IP  methods  (see  [68])  we  set  ADG(-,  ■)  =  A(-,  •)  where: 

A{uh,w)  =  (/eViih,  Vw)rh  -  [w])£h  +  0([it/J,  i 

+  (aeh~1Keluh},lw})£h,  V uh,  w  g  VDG  . 

The  weight  (3  =  {/3e}eg£-o,  depends  on  the  coefficient  k  and  therefore  it  might  vary  over  all 
interior  edges/faces.  For  any  e  G  £%  with  e  =  dT+  n  dT~ ,  we  set  /5e  as  follows: 


fie 


K 

K+  +  K~ 


where  =  «|  ± , 


and  we  define  the  coefficient  ne  as  the  harmonic  mean  of  k+  and  k  : 


2k+k 
+  KW 


(2.6) 


(2.7) 


Remark  2.1.  We  remark  that  one  could  take  Ke  as  min{K+,K  }  since  both  are  equivalent: 


min  {k+,  k  }  <  Ke  = 


2k+ 


K 


K+  +  K 


—  <  2min{K+,K  }  <  . 


(2.8) 


The  equivalence  relations  in  (2.8)  show  that  the  results  on  spectral  equivalence  and  uniform 
preconditioning  given  later  on  for  (2.5)  with  ne  defined  in  (2.7)  (the  harmonic  mean),  will 
automatically  hold  for  method  (2.5)  with  Ke  :=  min  {ft+,  /c-}.  To  fix  the  notation  and  to 
simplify  the  presentation  we  use  the  harmonic  mean  as  ne  [see  (2.7)]. 


The  symmetric  method  was  first  considered  in  [68]  and  later  in  [37,  Section  4]  for  variable 
coefficient  (although  there  it  was  written  using  a  slightly  different  notation  and  DG  was  only 
used  in  the  skeleton  of  the  partition).  It  was  later  extended  to  advection-diffusion  problems  in 
[25]  and  [32],  In  (2.5),  0  =  —  1  gives  the  SIPG(/3);  6  =  1  leads  to  NIPG(/3);  and  6  =  0  gives  the 
IIPG(/3)  discretizations.  The  penalty  parameter  ae  >  0  is  set  to  a  positive  constant;  and  for 


6 


B.  AYUSO  DE  DIOS,  M.  HOLST,  Y.  ZHU,  AND  L.  ZIKATANOV 


0^1  should  be  taken  large  enough  to  ensure  coercivity  of  the  corresponding  bilinear  forms.  We 
also  introduce  the  corresponding  family  of  IP (/3)-0  methods,  which  use  mid  point  quadrature 
rule  for  computing  all  the  integrals  in  (2.5).  That  is,  we  set  ADG(-,  •)  =  -4o(’>  •)  with 

A0(uh,w)  =  (a iVuh,Vw)Th  -  ({AeVti/Ja,,  [w])£h  +«(Iu],{kV%)£|1 
+  (aeh~1Ke'P^(luhj),V^(lwj))£h,  Vuh,  w  <E  VDG  . 

where  V®  :  L2{Sh)  P°(5/l)  is  the  L2-projection  onto  the  piecewise  constants  on  £/,.  We  note 
that  this  projection  satisfies  \\V® \\L2(£h)  =  1- 

Weighted  Residual  Formulation.  Following  [21]  we  can  rewrite  the  two  families  of  IP  meth¬ 
ods  in  the  weighted  residual  framework:  for  all  Uh,  w  6  VDG , 


A(uh,w)  =  (— V  •  (. KVuh),w)Th  +  (lKVuh],lw}i-pe)£o  +  (luh],Bi(w)}£h,  (2.10) 

A0(uh,w)  =  (-V  •  ( KVuh),w)Th  +  +  {luhj,V^(B1(w)))£h,  (2.11) 

where  B\  is  defined  as: 

Bi(w)  =  OJ[KVwJpe  +  aeh~1Kelw},  V e  e  £h.  (2.12) 

Throughout  the  paper  both  the  weighted  residual  formulation  (2.10)-(2.11)  and  the  standard 
one  (2.5)-(2.9)  will  be  used  interchangeably. 


Continuity  and  Coercivity.  The  family  of  methods  (2.5)  and  (2.9)  can  be  shown  to  provide 
an  accurate  and  robust  approximation  to  the  solution  of  (1.1).  We  define  the  energy  norm 

IIMIIdg1 

IIKIIllo  :=  E  V*(/illo,T  +  E  kA-'IIMo,.-  <2-13) 

T&Th  e&£h 

For  the  classical  IP  methods,  the  bilinear  form  «4(-,-)  is  continuous  and  coercive  in  the  above 
norm,  with  constants  independent  of  the  mesh  size  h  and  the  coefficient  k: 

Continuity:  \A{uh,w)\  <  |||u^|||DG  ||MIIdg>  V uh  ,  w  G  VDG  (2.14) 

Coercivity:  A(uh,uh)  >  \\\uh\\\2DG  ,  \/uh  £  VDG  .  (2.15) 

The  proof  of  (2.15)  and  (2.14)  is  standard  and  could  be  found  in  [37].  We  sketch  it  here  for 
completeness.  Note  first  that  for  each  e  £  £fp  such  that  e  =  dT+  n  dT~ ,  the  weighted  average 
can  be  rewritten  as: 

{nVu}0e  =  /3e(n+(Vu)+)  +  (1  -  /3e)(/W(V«)-) 

=  _^+(Vu)+  +  h  _/c~(Vm)~ 

K,'  T  Av  K'  T  Av 

=  [(Vu)+  +  (V«n  =  Ke{{Vu}}  .  (2.16) 

The  trace  inequality  [1],  the  inverse  inequality  [28],  and  (2.8)  then  imply  the  following  bounds 


Bell'S  ftVw,]}  pe  Ho,e  —  Ct(ne)2  (^ll^'ullo,T+UT-  +  h2  |Vu| \^+YjT~  ) 

<  2(ACe)C't(l  +  0  (ac+||Vu||2t+  +«-||V«||g)T-) 


(2.17) 
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Then,  using  Cauchy-Schwarz  inequality,  the  estimate  (2.17),  and  the  arithmetic-geometric  in¬ 
equality  we  have 


Kf/eVuJ^,  H)sJ 


KefVu|[u;]c?s 


< 


8Ct(l  +  CgJ 

Qte 


Y  Kt||Vu||o,T+  \  ^^eVlIMIIo.e 

T&Th  e&£h 


Then,  (2.14)  follows  from  Cauchy-Schwarz  inequality  and  the  estimate  (2.18).  The  inequality 
(2.15)  is  proved  by  setting  w  =  u  in  (2.5)  and  taking  into  account  the  estimate  (2.18)  (also  with 
w  =  u).  We  have, 


A(u,u)=  Y  KT\\Vu\\\T  +  aeY  «e^e1HNIIo,e-(l-^X-|[«V«Jy9e,[«])£h 
TgTh  e&£h 

>  IIMII DG  ~  I1  “  °\  l({{KV^/3e>  N>fJ 


>  I  1  -  8C‘(1  +  C,L)1  E  «rl|V«||g,r  + 


ae 


TeTh 


cx.e  '  nehe  | 
ee£h 


U 


l0,e  > 


and  (2.15)  follows  immediately  by  taking  ae  large  enough  (if  9  A  1).  Moreover,  notice  that 
both  constants  in  (2.14)  and  (2.15)  depend  on  the  shape  regularity  of  the  mesh  partition  but 
are  independent  of  the  coefficient  k. 

For  the  IP-0  methods  (2.9),  similar  properties  can  also  be  easily  shown  to  hold  following  the 
same  arguments,  albeit  in  a  different  energy  norm,  defined  as: 


u 


\dG0  KtI|Vu||o)T  +  Y,  Kehe  1||T’e  (M)||o,e- 

TgTh  e&£h 


(2.18) 


For  both  families  of  methods,  optimal  error  estimates  in  the  energy  norms  (2.13)  and  (2.18)  can 
be  shown,  arguing  as  in  [5].  See  also  [8]  for  further  discussion  on  the  L2-error  analysis  of  these 
methods. 

We  now  establish  the  spectral  equivalence  between  M(-,  •)  and  Mo(-,  •)• 


Lemma  2.2.  Let  A(-,-)  be  a  bilinear  form  corresponding  to  a  Type-1  IP  (ft)  method  (2.5)  and 
let  Aq ( • ,  ■ )  be  the  corresponding  Type-0  bilinear  form  as  defined  in  (2.9).  Then  there  exists  a 
positive  constant  cq  =  co(a),  depending  only  on  the  shape  regularity  of  the  mesh  and  the  penalty 
parameter  a  (but  independent  of  the  coefficient  n  and  the  mesh  size  h )  such  that, 

Ao(u,u)  <  A(u,u)  <  co(a)Ao(u,u)  \/u  G  VDG.  (2-19) 

Proof.  The  lower  bound  follows  immediately  from  the  fact  that  the  projection  V '*?  is  an  L2(e)- 
orthogonal  projection  and  therefore  has  norm  1. 

To  show  the  upper  bound  it  is  enough  to  show  that 

Y  “e^rVlIMIIo.e  «rl|Vu||o,r  +  aeKlKe\\VeM\\l,e)-  (2-20) 

Te.'Tk 

Adding  and  subtracting  V®  [u]  in  the  term  on  the  left  side  above  and  using  again  that  V °  is  the 
L2(e)-orthogonal  projection  on  the  constant  functions  we  have  for  each  face  e, 

II Hilo, e  =  IT,°(H)llh  +  II M  -  Pe°(H)ll5,e  • 


(2.21) 
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Hence,  we  only  need  to  estimate  the  last  term  on  the  right  side  of  the  above  inequality  and  show 
that  this  term  is  bounded  by  the  right  side  in  the  inequality  given  in  (2.20).  Observe  that  on 
each  e  G  £,  e  =  8T+  0  dT~ , 

II H  -  A°(M)IIS,,  <  ll«+  -  AVOIft,  +  ll«-  -  A°(<OllL . 

In  [9,  Lemma  2.1]  it  was  further  shown  that 

KX\\U±  -^e(u±)\\o,e  £  \u±\l, T±  •  Ve  e  £h 

Thus,  multiplying  the  above  expression  by  the  harmonic  average  Ke  and  taking  into  account  (2.8) 
we  arrive  at 

Y  “e/»r1«e||N||o,e  <  C(a )  Y  Kt\Nu\\o,T+  Y  V||Te°H  llo,e  > 
eESh  Te  Th  e€£h 

which  concludes  the  proof.  □ 

3.  Space  decomposition  of  the  Vdg  space 

In  this  section  we  introduce  a  decomposition  of  the  -space  that  will  play  a  key  role  in  the 
design  of  the  solvers  for  the  DG  discretizations  (2.5)  and  (2.9).  In  [9]  (see  also  [24])  it  is  shown 
that  the  discontinuous  piecewise  linear  finite  element  space  VDG  admits  the  decomposition: 
VDG  =  VGR  ©  Z  where  VGR  denotes  the  standard  Crouziex-Raviart  space  defined  as 

VGR  =  {«£  L2(H)  :  v\T  G  P  1{T)  VT  G  Th  and  Te0([u]  •  n)  =  0  Ve  €  £°h]  ,  (3.1) 

and  the  complementary  space  Z  is  a  space  containing  functions  having  zero  averages  at  the 
midpoints  of  the  internal  edges: 

Z={z€  L2(n)  :  Z\T  G  P \T)  VT  G  Th  and  V°fv}  =  0,  Ve  G  ££}  . 

In  [9]  it  was  shown  that  this  decomposition  satisfies  Aq(v,z)  =  0,  for  all  v  G  VGR  and  z  G  Z. 
We  now  modify  the  definition  of  Z  above  in  order  to  account  for  the  presence  of  a  coefficient  in 
the  problem  (1.1);  we  define 

Zp  =  {z  G  L2(H)  :  zjT  G  Px(T)  VT  G  Th  and  V°e({{z} =  0,  Ve  G  £°h}  ,  (3.2) 

where  the  weight  ft  was  defined  earlier  in  (2.6).  Note  that  the  weight  (5  depends  on  the  coefficient 
k,  and,  as  a  consequence,  the  space  Zp  is  also  coefficient  dependent.  In  what  follows,  we  shall 
show  that  Zp  is  a  space  complementary  to  VGR  in  VDG  and  the  corresponding  decomposition 
has  properties  analogous  to  the  properties  of  the  decomposition  VDG  =  VGR  ©  Z  given  in  [9] 
for  the  Poisson  problem. 

For  any  eG4  with  e  C  T  G  Th,  let  <pe,T  be  the  canonical  Crouzeix-Raviart  basis  function  on 
T,  which  is  defined  by 

Te,T|r  e  P1(T),  <pe,T(me /)  =  4,e'  Ve'  G  £h(T),  and  <pe,T(x)  =  0  Vx  0  T, 

where  me  is  the  mass  center  e.  We  will  denote  ut  and  ng  as  the  number  of  simplices  and  faces 
(or  edges  when  d  =  2)  respectively.  We  also  denote  ube  as  the  number  of  boundary  faces. 

Proposition  3.1.  For  any  u  G  VDG  there  exists  a  unique  v  G  VGR  and  a  unique  zpe  G  Zp  such 
that  u  =  v  +  zpe  ,  that  is 

VDG  =  V£R  ©  Zp.  (3.3) 

Proof.  Throughout  the  proof,  for  simplicity,  let  us  set  /3+  =  j3e,  /3~  =  (1  —  /3e),  and  ipf  =  ipeT± 
for  any  eG^  with  e  =  dT+  n  dT~ .  We  also  denote  =  <^e,T  for  any  e  G  £ft  with  e  =  8T  n  d£l. 
Note  that 

dim  VDG  =  (d  +  1  )iit  =  2  ue  —  ube, 


MG  FOR  DG  WITH  JUMP 


9 


and  it  is  also  obvious  that  {cp±}ee£o  U  {fe}e&£a  form  a  basis  for  VDG .  Notice  that  (3+  +  (3  =  1, 
therefore  we  can  express  any  u  £  VDG  as 


u(x )  =  ^  U+ (m,e)yt  (.x)  +  ^2  U  (me)Pe  ( x )  +  u(me)<pe(x) 

eef°  ee£°  ee£S 

=  ^{P~u+{me)  + P+u-(me))(vt{x)  +  Ve(x)) 

+  ^2(u+(me)  -  u~(m,e))(/3+ip+(x)  -  f3~ip-(x))  +  ^2  u(me)<pe(x) 

e££h  ee£f 

=  (u  [ (<pt (x)  +  <p7 (x)) 

e&£°  I  Je  ' 

+  fn  / Nn+ds)  03+¥>e  (*)  -  P~Ve{x))  +^2  (t~\  f  MMS)  ^e(a 

ocfo  Vlel  Je  /  Vlel  Je  / 


°g£° 

=  v{x)  +  zp(x). 


26 £ ? 


Then  for  each  e  €  ££,  we  set 

:=  ^ (*)  +  (3-4) 

and  let 

VM  ■■=  P+vi(*)  -  P'vztx)  =  {  .  (3-5) 

and  ijje(x)  :=  0  for  all  x  0  T+  UT_.  In  the  definition  (3.5)  of  ip^(x),  we  have  used  <^“(x)  =  0  for 
x  £  T+  and  <p£(x)  =  0  for  x  £  T~ .  Finally,  when  e  £  £fr  with  e  =  dT  fl  d£l  for  some  T,  we  set 

i>l(x)  =  <pe(x),  Vx  £  T.  (3.6) 


It  is  then  straightforward  to  check  that 

VhR  =  span{ipGR}ee£o,  and  Zp  =  span{^}ee^. 
Hence,  for  all  u  £  VDG  there  are  unique  v  £  VGR  and  zpe  £  Zg  defined  by 


such  that  u  =  v  +  zg.  This  shows  (3. 3) and  concludes  the  proof.  □ 


Remark  3.2.  As  we  pointed  out  in  the  introduction,  the  definition  of  the  subspace  Zg  clearly 
depends  on  the  coefficient  k  since  /3  depends  k.  Such  dependence  is  also  often  seen  in  algebraic 
multigrid  analysis,  where  the  coarse  spaces  depend  on  the  operator,  and  are  in  fact  constructed 
in  this  way,  the  aim  being  to  increase  robustness  of  the  methods. 

In  the  proof  of  Proposition  3.1  above  we  have  introduced  the  basis  in  both  VfGri  and  Zg.  The 
canonical  Crouzeix-Raviart  basis  functions  are  denoted  with  {<pGR }ee£°  f°r  non-conforming  and 
these  functions  are  continuous  at  the  mass  centers  me  of  the  faces  e  £  £ The  basis  in  Zg, 
{ipeJe^Sh  consists  of  piecewise  P1  functions,  which  are  discontinuous  across  the  faces  in  £}x.  In 
fact,  for  any  z  £  Zg  such  that  z  =  Ylee£h  z^e  with  ~e  £  1R  for  any  e  £  £g,  we  have 

(|z]n+)(me/)  =  ze>,  Ve'  £  £h. 
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To  see  this,  evaluating  the  jump  of  z  at  mei  gives 

(Mn+)("v)  =  Ze(Hein+)(me')  =  Ze'(bPe'ln+)(me>) 

e€£h 

(  zet (/3ei  -  (fie,  -  1))  =  ze>,  e'  G  £%, 

1  Ze'i  e'  G  £%. 

This  relation  will  also  be  used  later  to  obtain  uniform  diagonal  preconditioners  for  the  restric¬ 
tions  of  A(-,  •)  and  Ao(-,  •)  on  Zp. 

Remark  3.3.  For  mixed  boundary  value  problems,  that  is,  8£l  contains  both  Neumann  boundary 
Tjv  0  and  Dirichlet  boundary  T jj  with  d£l  =  T/j  U  T^r,  the  definition  of  the  basis  functions  on 
the  boundary  faces  [see  (3.6)]  need  to  be  changed  as: 

ipeR(x)  =  ipe,T{x),  e  =  dTnTN,  for  all  x  G  T,  ,  . 

ipe(x)  =  <pe,T(x),  e  =  dTnTn,  for  all  iGT,  ' 

Thus,  in  case  T^r  0  the  dimension  of  VGR  is  increased  (by  adding  to  it  functions  that 
correspond  to  degrees  of  freedom  on  Tat)  and  the  dimension  of  Zp  is  decreased  accordingly. 
Clearly  things  balance  out  correctly:  the  identity  VDG  =  VGR®Zp  holds,  and  also  the  analysis 
carries  out  with  very  little  modification. 

Next  lemma  is  a  simple  but  a  key  observation  used  in  the  design  of  efficient  solvers  and 
preconditioners  and  it  shows  that  the  restriction  of  Aq(-,  ■)  :  VDG  x  VDG  i-A  R  to  VGR  x  Zp 
vanishes. 

Lemma  3.4.  Let  u  G  VDG  be  such  that  u  =  v  +  z  with  v  G  VGR  and  z  G  Zp.  Let  »4o(-)  •)  be 
the  bilinear  form  defined  in  (2.9).  Then, 

A0{v,z)  =  0  VuGT4cr,  M z  &  Zp.  (3.8) 

Furthermore  if  Aq(-,  •)  is  symmetric,  then  Aq{v,z)  =  Aq(z,v)  =  0  for  all  v  G  VGR,  and  for  all 
z  G  Zp  and  the  decomposition  (3.3)  is  Aq- orthogonal,  namely,  VGR  T^0  Zp. 

Proof.  From  the  weighted-residual  form  of  „4o(v)  given  in  (2.11),  for  all  v  G  VffR,  and  all 
z  G  Zp  we  easily  obtain 

Afiv,z)  =  (-V  •  {nVv),z)Th  +  <[KV«Mzh_A>£o  +  M^Biiz))) £h  =  0. 

In  the  equation  above,  the  first  term  is  zero  due  to  the  fact  that  v  is  linear  in  each  T  and  that 
the  coefficient  k  is  a  constant  on  T;  the  second  term  vanishes  from  the  definition  of  Zp  (since 
[kVu]  is  constant  on  each  e  G  £f)  and  last  term  vanishes  as  well  independently  of  the  choice  of 
6  (or  equivalently  the  choice  of  B\{v))  from  the  definition  of  the  space  VPR.  Moreover,  in  the 
case  when  Mo(-,  •)  is  symmetric  and  positive  definite  we  have  that  Ao(v,z)  =  Ao{z,v),  for  all 
v  G  VGR  and  for  all  z  G  Zp.  Thus,  for  the  symmetric  method  Mo(-,  •),  the  spaces  VGR  and  Zp 
are  indeed  Mo-orthogonal.  The  proof  is  complete.  □ 

4.  Solvers  for  IP-0  methods 

In  this  section  we  show  how  Proposition  3.1  and  Lemma  3.4  can  be  used  in  the  design  and 
analysis  of  uniformly  convergent  iterative  methods  for  the  IP-0  methods.  We  follow  the  ideas 
and  analysis  introduced  in  [9]  and  point  out  the  differences.  We  first  consider  the  approximation 
to  problem  (1.1)  with  ADG(-,  •)  =  Aq(-,  •)•  To  begin,  let  Aq  be  the  discrete  operator  defined  by 
( Aqu,w )  =  Ao{u,w)  and  let  Ao  be  its  matrix  representation  in  the  new  basis  (3.7).  We  denote 
by  u  =  [z,v]T,  f  =  [fz,  fv]7  be  the  vector  representation  of  the  unknown  function  u  and  of  the 
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right  hand  side  /,  respectively,  in  this  new  basis.  A  simple  consequence  of  Lemma  3.4  is  that 
the  matrix  Aq  (in  this  basis)  has  block  lower  triangular  structure: 


Aq  = 


A;f  0 
Ag*  Ag” 


(4.1) 


where  Aq2,Aqu  are  the  matrix  representation  of  Aq  restricted  to  the  subspaces  Zg  and  VCR, 
respectively,  and  Aq2  is  the  matrix  representation  of  the  term  that  accounts  for  the  coupling 
(or  non-symmetry)  Ao(ipz,  <pCR).  As  remarked  earlier,  for  Type-0  SIPG(/3)  discretization,  the 
stiffness  matrix  Ao  is  block-diagonal. 

Here,  we  have  given  a  2D  example,  with  two  squares  Hi  =  [— 0.5,0]2  and  D2  =  [0,0. 5]2  inside 
the  domain  D  =  [—1,  l]2.  We  set  the  coefficients  k(x)  =  1  for  all  1  G  Sli  U  SI2  and  n(x)  =  10~3 
when  x  G  D  \  (Hi  U  f^),  see  Figure  4.1.  We  choose  the  penalty  parameter  ae  =  20  in  (2.9). 
Figure  4.2  shows  the  sparsity  patterns  of  the  Type-0  IP(/3)  methods  with  standard  nodal  basis, 


G2 

Q1 

Figure  4.1.  Computational  domain  and  unstructured  mesh. 


while  Figure  4.3  shows  the  sparsity  patterns  of  the  stiffness  matrices  for  the  IP-0  (/3)  methods 
after  changing  from  the  standard  nodal  basis  to  the  basis  (3.4)-(3.5)  induced  by  the  splitting 

% 


VDG =  yCR 


SIPG  (nnz=  10834) 


NIPG  (nnz=1 0834) 


IIPG  (nnz  =  8886) 


Figure  4.2.  Non-zero  pattern  of  the  matrix  representation  in  the  standard 
nodal  basis  of  the  operators  associated  with  Type-0  IP  methods.  From  left 
to  right:  SIPG,  NIPG  and  IIPG  methods. 

Clearly,  as  in  the  constant  coefficient  case,  a  simple  algorithm  based  on  a  block  version 
of  forward  substitution  provides  an  exact  solver  for  the  solution  of  the  linear  systems  with 
coefficient  matrix  Ao-  A  formal  description  of  this  block  forward  substitution  is  given  in  the 
next  Algorithm. 
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Figure  4.3.  Non-zero  pattern  of  the  matrix  representations  according  to  the 
basis  (3.3)  of  the  operator  associated  with  Type-0  IP  methods;i.e.,  Ao-  From  left 
to  right:  SIPG,  NIPG  and  IIPG  methods. 


Algorithm  4.1.  1.  Find  z  G  Zp  such  that  Ao(z,ip)  =  (f,ip)Th  f°r  aH  ^  £  Zp 

2.  Find  v  G  VCR  such  that  Aq(v,  ip)  =  (/,  p)rh  —  Ap{z,  p)  for  all  P  £  VCR 

3.  Set  u  =  z  +  v 

Notice  that  the  above  algorithm  requires  the  solution  of  Aq(-,  •)  on  Zp  (Step  1.  of  the  algo¬ 
rithm)  and  the  solution  of  Ao(-,-)  on  V/fR  (Step  2.  of  Algorithm  4.1).  Unlike  the  situation  in 
[9],  due  to  the  jump  coefficient  in  (1.1),  the  solution  on  VCR  is  more  involved,  and  therefore  we 
postpone  its  discussion  and  analysis  until  Section  5.  On  the  other  hand,  as  we  show  in  the  next 
section,  the  solution  on  Zp  (the  space  regarded  as  containing  “high  frequency”  components) 
can  be  handled  efficiently  using  CG  with  a  diagonal  preconditioner. 

4.1.  Solution  on  Zp .  The  first  result  in  this  section  establishes  the  symmetry  of  the  restrictions 
of  the  bilinear  forms  (of  both  Type-0  and  Type-1)  on  Zp. 

Lemma  4.2.  Let  A(-,-)  be  the  bilinear  form  of  a  non-symmetric  Type-1  IP  method  as  defined 
in  (2.5)  and  let  Ao(-,-)  be  the  corresponding  Type-0  bilinear  form.  Then  the  restrictions  to  Zp 
of  both  Aq{-,  •)  and  A(-,  •)  are  symmetric.  Namely,  for  9  =  —1,  0,  1  and  for  all  z  G  Zp  and 
4>  G  Zp  we  have 

Aq(z,  <f>)  =  Ao((f,  z)  and  A{z,<f)  =  A{(j>,  z). 

Proof.  If  9  =  —  1  there  is  nothing  to  prove,  since  in  this  case  both  bilinear  forms  are  symmetric. 
Hence  we  only  consider  the  cases  9  =  0  or  9  =  —1.  Integrating  by  parts  and  using  the  fact  that 
z  G  Zp  and  6  Zp  are  linear  on  each  element  T  shows  that 

0  =  (-V  ■  («W),  Vz)Th  =  {nV^,Vz)Th  -  ({{^VV^,  \z\)Sh  ~  ([«V^],  Mi-zUfj’  ■  (4-2) 

Hence,  from  the  definition  (3.2)  of  the  Zp  space,  it  follows  that 

(k,VV>,  Vz)Th  =  {{n\/ip}pe,  lzj)£h  =  {{K\7z}pe,  m)sh  yz,tpe  Zp.  (4.3) 

Substituting  the  above  identity  in  the  definition  of  the  bilinear  form  (2.9)  then  leads  to 

A0(z,ip)  =  0(H,1kV%)4  +  (V°(lz}),Kem)£h 

=  9(KVlf,Vz)Th  +  ('Pe(m),Kelzf)£h  =  A0(lf,z). 

This  shows  the  symmetry  of  Ao(-,  •)  on  Zp.  The  symmetry  for  A(-,  •)  on  Zp  then  follows  from  the 
result  for  Ao(-,  ■),  since  the  difference  between  these  bilinear  forms  is  obviously  symmetric.  □ 
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We  now  study  the  conditioning  of  the  bilinear  forms  A(-,  •)  and  *4o(u  •)  on  Zp.  For  all  z  G  Zp, 
and  for  all  f>  G  Zp  with 

z  =  E  zeil>ze  G  Zp,  and  f  =  E  fe ip*  G  Zp. 

we  introduce  a  weighted  scalar  product  (•,  •)*  :  Zp  X  Zp  1R  and  the  corresponding  norm  ||  •  ||*, 
defined  as  follows 


(z,<P)*  ■=  E  ]KeZ‘ 


e&£h 


z  *  :=  [z,z 


(4.4) 


Observe  that  the  matrix  representation  of  the  above  weighted  scalar  product  (in  the  basis  given 
in  (3.5)),  is  in  fact  a  diagonal  matrix.  The  next  result  shows  that  the  restrictions  of  «4(-,  •)  and 
A(\(-7  •)  to  Zp  are  spectrally  equivalent  to  the  weighted  scalar  product  (•,•)*  and  therefore  their 
matrix  representations  are  spectrally  equivalent  to  a  diagonal  matrix. 

Lemma  4.3.  Let  Zp  be  the  space  defined  in  (3.2).  Then  for  all  z  G  Zp,  the  following  estimates 
hold 

\\z\\l  <  A0(z,  z)  <  \\z\\l,  (4.5) 

and  also 


A{z,  z) 


(4.6) 


Proof.  Let  us  fix  z  G  Zp,  z  =  Yle££h  Zef’e-  From  the  definition  of  V® (H)||o,e>  ^  immediate 
to  see  that 


l0,e 


=  ez„ 


>0'W)IIS,« 


E 


Kp  —  Z^. 


h, 


(4.7) 


Thus,  we  have  that 

E  Kehei§'Pe 

e&£h  ee£h 

From  this  relation  it  is  easy  to  show  both  estimates  (4.5)  and  (4.6).  First,  to  show  (4.5),  we 
notice  that  taking  into  account  (4.3)  together  with  the  estimate  (2.8),  it  follows  that 

E  «tI|V*||§,t  =  (kV*,  Vz)Th  =  ({{^z}pe,  {z\)£h  =  (KelVz},V°e(lzj))£h 


TeTh 


1/2 


< 

r-j 


E  Kt\\Vz\\It  E  Ke\\he1/2'Pe(lz. 


,T€Th 


e(z:Eh 


and  therefore, 


I2  =  IMI2 

I0,e  11^11*5 


(4.8) 


E  ^llVzHo.T  £  E  KeWke  1/2pe(W 

Te77  e&£h 

and  since  z  G  Zp  was  arbitrary,  we  have  that  Ao(z,z)  <  ||^||*,  for  all  z  G  Zp  and  this  proves 
the  upper  bound  in  (4.5). 

To  prove  the  lower  bound,  we  use  the  coercivity  estimate  for  the  bilinear  form  Aq ( • ,  •)  in  the 
energy  norm  ||| •  ||lz?G0  [see  (2-18)].  For  all  z  G  Zp  we  have 


A0(z,z) 


> 


Idgo  —  E  +  E  MV1/2P?(M)llo,e 

TtzTh  e^Eh 

|2 


>  E  ^\\Kll2r0e(l4)\\l,e  =  114,- 

e€.£h 

Since  this  is  the  desired  bound  we  conclude  the  proof  of  (4.5). 

The  proof  of  (4.6)  follows  easily  from  the  spectral  equivalence  between  A(-,  ■)  and  A(j(-,  ■) 
given  in  Lemma  2.2.  □ 
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Observe  that  last  result  guarantees  that  the  linear  systems  on  Zg  can  be  efficiently  solved 
by  preconditioned  CG  (PCG)  with  a  diagonal  preconditioner.  As  a  corollary  of  the  result  in 
Lemma  4.3,  the  number  of  PCG  iterations  will  be  independent  of  both  the  mesh  size  and  the 
variations  in  the  PDE  coefficient. 

We  now  show  that  in  the  particular  case  of  IIPG-0  method,  the  matrix  representation  of 
*4o('>  •)  on  -Z/j  is  itself  a  diagonal  matrix. 

Lemma  4.4.  Let  Ao(-,  •)  be  the  bilinear  form  of  the  non-symmetric  IIPG{f3 )  Type-0  method 
(2.9)  with  9  =  0.  Let  {'ifz}e^£h  be  the  basis  for  the  space  Zp  as  defined  in  (3.5).  Let  Aqz  be  the 
matrix  representation  in  this  basis  of  the  restriction  to  the  subspace  Zp  of  the  operator  associated 
to  -4o(’i  •)•  Then,  kfz  is  diagonal. 


Proof.  Note  that  from  the  definition  (2.9)  of  the  method  (9  =  0)  together  with  (4.3)  we  have 

MzA)  =  (Vz,  WOr,  -  <«VzU,  M)£h  +  (aeK^eV^lzD^dmSn 

=  (aeK1KeV®(lz}),V°(l'if}))£h  \/z,ifeZp  .  (4.9) 

Let  {V;f}{ee£?l}  be  the  basis  functions  (3.5).  To  prove  that  Aqz  is  diagonal  it  is  enough  to  show 
that  for  the  basis  functions  (3.5),  the  following  relation  holds: 

AoiVe,  V’eO  =  Ce<5e,e'  Ce  /  0,  V  e  €  £h  ,  (4.10) 

where  5e>e/  is  the  delta  function  associated  with  the  edge/face  e.  We  now  show  (4.10).  Observe 
that  the  supports  of  ip*  and  ip*,  have  empty  intersection  unless  e,  e'  C  T  for  some  T  £  Th-  Let 
T  n  dQ  =  0  be  an  interior  element,  then  from  (4.9)  and  the  mid-point  integration  rule,  we  have 


•AoiVeiVe’)  =  otehel 


aehe  1Ke[2ip*(me)][2ipf,(me)\ 


=  Aaehe  1ne5e^ei  e,  e'  C  dT  e,  e'  £  , 

which  shows  (4.10)  for  interior  edges  with  ce  =  4ae/i“1Ke.  For  boundary  edges/faces  the  con¬ 
siderations  are  essentially  the  same  and  therefore  omitted.  The  proof  is  complete  since  the 
relation  (4.10)  readily  implies  that  the  off-diagonal  terms  of  Aq2  are  zero.  □ 


We  remark  that  this  lemma  will  be  essential  for  proving  the  uniform  convergence  of  the 
iterative  solvers  for  the  non-symmetric  methods. 


5.  Robust  preconditioners  on  VffR 

In  this  section  we  develop  efficient  and  robust  solvers  for  the  solution  in  the  CR  space.  By 
virtue  of  the  spectral  equivalence  between  A  and  „4o  given  in  Lemma  2.2,  it  is  enough  to  focus 
on  the  development  of  solvers  for  „4o  restricted  to  the  subspace  and  this  is  what  we  do  next. 
Observe  that  the  restriction  of  A<j  ( • ,  • )  to  the  subspace  VCR  reduces  to  the  P1-nonconforming 
finite  element  discretization  of  (1.1) 

Find  u  £  :  Ao(u,w)  =  (Wu,  Vw)rft  =  (/,  w)  \/ w  £  VffR.  (5.1) 

We  denote  Aq  r  as  the  operator  induced  by  (5.1).  For  the  analysis  that  follows,  we  need  the 
following  semi- norms  and  norms  for  any  v  £  VPR: 

Ml ,h,K  :=  KTl|V'y||o,T  j  Ml, ft, Hi  :=  Iloilo, T  >;■  (5-2) 

TeTh  T&Th  ,  Tcat 

M 

II  ll 2  I  II  1 1 2  II  1 1 2  II  1 1 2  ,  |  1 2  (r-  q \ 

ItIIoa  •  -  2^KlnillvlloA  »  \\v\\i,h,K  Mllo,K  +  \v\i,h,K  •  (5-3) 

i=  1 

If  the  coefficient  k(x)  is  a  constant  or  a  smooth  function;  one  can  find  some  uniform  precon¬ 
ditioners  in  the  literature  (see  the  references  in  [9]).  The  case  of  jump  coefficients  k  has  also 
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been  considered  by  several  authors  [65],  [64],  with  the  assumption  on  the  quasi- uniformity  in 
the  coefficients. 

Since  (5.1)  is  a  symmetric  problem,  from  the  classical  theory  of  PCG  we  know  that  the  con¬ 
vergence  rates  of  the  iterative  method  with  preconditioner,  say  B,  for  Aq  are  fully  determined, 
in  the  worst  case  scenario,  by  the  condition  number  of  the  preconditioned  system:  K(BAqR). 
However,  for  problems  with  large  jumps  in  the  coefficient  k,  it  has  also  been  observed  in  [49,  78] 
that  the  spectrum  o(BAq  r)  might  contain  a  few  very  small  eigenvalues,  which  result  in  a  ex¬ 
tremely  large  value  of  K(BAqR )  (independently  of  how  fine  is  the  mesh),  but  they  seem  to  have 
no  influence  on  the  efficiency  of  the  preconditioner  and  the  overall  convergence  of  the  iterative 
method.  Therefore,  following  the  approach  introduced  in  [78],  we  consider  the  asymptotic  con¬ 
vergence  rate  of  the  PCG  algorithm,  which  is  determined  by  the  so-called  effective  condition 
number,  as  defined  below  in  Definition  5.1. 

The  setting  that  we  need  for  the  analysis  require  some  more  preliminaries.  Let  us  first  briefly 
review  few  standard  results  for  the  PCG  algorithm.  Suppose  that  the  spectrum  of  BA//r, 
o(BAqR),  is  divided  in  two  sets:  (j(BAqR )  =  MBACX)  U  o-i(BAqR),  where  MBA°R)  = 
{Ai,...,Amo}  contains  of  all  very  small  (often  referred  to  as  “bad”)  eigenvalues  and  the  re¬ 
maining  eigenvalues  (bounded  above  and  below)  are  in  cf\(BAq  r)  =  {Amo+i, . . . ,  Ajv}  ;  that  is 
A j  €  [a,  b]  for  j  =  mo  +  1, . . . ,  N,  with  N  =  dim(VjpR)  =  ng  —  ube •  Then,  the  error  at  the  k- th 
iterate  of  the  PCG  algorithm  is  bounded  by  (see  e.g.  [6,  51,  7]): 


u-uk\\lXK<2(K(BACR)-l)m° 


( \fWa~ 1 
\  \/b/a  +  1 


k—m  o 


\U-  Uo\\lih,K 


(5.4) 


From  the  above  estimate,  one  may  conclude  the  following:  If  there  are  only  a  few  small 
eigenvalues  of  BAqR  in  (Tq(BAqR),  then  the  asymptotic  convergence  rate  of  the  resulting 


PCG  method  will  be  dominated  by  the  factor  ^=— i.e.  by  b/a  where  b  =  Xn(BAqR )  and 

a  =  Am0+i(-£L4Q  R)-  The  quantity  (b/a)  which  determines  the  asymptotic  convergence  rate  is 
often  called  effective  condition  number  and  its  formal  definition  follows. 


Definition  5.1.  Let  V  be  a  real  IV-dimensional  Hilbert  space,  and  21  :  V  — >  V  be  a  symmetric 
positive  definite  linear  operator,  with  eigenvalues  0  <  Ai  <  •  •  •  <  Xjy.  The  vriQ-th  effective 
condition  number  of  21  is  defined  by 


ICmo+ i(2l) 


Ay  (21) 
Amo+l  (21) 


We  note  that  this  characteristic  was  previously  used  in  the  context  of  deflation  based  precon¬ 
ditioner  (for  determining  the  asymptotic  convergence  rate  of  PCG  for  semidefinite  systems  [45]). 

In  our  analysis  of  two-level  and  multilevel  preconditioners  for  problem  (5.1),  we  will  mainly 
focus  on  estimating  the  effective  condition  number.  Related  works  which  derive  estimates  on 
the  effective  condition  number  in  the  context  of  multigrid  and  domain  decomposition  methods 
for  conforming  discretizations  are  [49,  78,  80].  To  improve  the  flow  of  the  discussion  when 
presenting  the  analysis  of  the  multilevel  preconditioner  for  the  nonconforming  discretizations, 
we  start  by  introducing  a  two  level  method  and  provide  a  detailed  analysis  and  bounds  of  the 
condition  number  of  the  preconditioned  system. 


5.1.  Two  level  preconditioner  for  Alo(-,  •)  on  V/fR.  In  this  subsection  we  will  construct  an 
additive  preconditioner  (parallel  subspace  correction  method)  with  a  conforming  space  V~oni  C 

V/fR  as  a  coarse  space  plus  a  pointwise  relaxation  (point- Jacobi  or  Gauss-Seidel  method).  More 
precisely  to  define  the  two  level  parallel  subspace  correction  method,  we  consider  the  following 
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overlapping  space  decomposition  of  V^R: 

VfR  =  VfR  +  H~conf.  (5.5) 

Here,  we  allow  for  different  h:  one  could  set  h  =  h,  or  choose  a  coarser  approximation  i.e., 

h  =  H  >  h.  In  addition  we  will  assume  here  that  the  meshes  77  and  7r  are  nested.  On  IAconf 

n  h 

we  consider  the  standard  conforming  P1-approximation  to  (1.1):  Find  x  £  V~ont  such  that 

a(x,  rj)  :=  A0(x,  v)  =  [  ■  Vr/cte  =  (/,  77),  V??  G  H~conf  .  (5.6) 

The  bilinear  form  given  above  defines  an  “energy”  scalar  product  and  its  naturally  induced 
weighted  semi-norm  is: 

\x\Ik,d  ■=  [  ^X\2dx  ,  \/XeH\D),  Den.  (5.7) 

J  D 

We  next  discuss  the  matrix  form  of  the  additive  preconditioner  and  then  its  operator  form. 
The  reason  for  discussing  both  these  equivalent  forms  is  that  the  matrix  form  is  used  to  im¬ 
plement  the  method,  while  the  operator  form  is  more  suitable  for  the  analysis  and  the  spectral 
equivalence  results  that  we  prove. 

5.1.1.  Definition  of  the  preconditioner:  Matrix  notation.  We  begin  by  introducing  the 
necessary  notation  pertaining  to  vectors  and  matrices.  We  write  Aq  r  =  Aqv  =  Do  —  D  —  L* 
where  Do?  and  L  are  the  diagonal  and  the  strict  lower  triangular  part  of  the  stiffness  matrix 
Aq  representing  the  restriction  of  Ho(-,  •)  on  R .  We  use  similar  notation  for  the  space  V~onf 
denoting  by  Ac  the  stiffness  matrix  representing  the  bilinear  form  a(-,  •)  on  the  conforming  space 
V~onf.  For  v  G  VfR,  v(x)  =  Yhe&s°  vetPeR(x )?  we  denote  by  v  G  RncH  the  vector  of  degrees  of 
freedom  of  v(x),  that  is  v  has  components  {ve}e&£°.  Here  ucr  =  ng  —  hbe  is  the  dimension 
of  the  space  VffR  and  we  shall  denote  by  nc  the  dimension  of  V~onf.  The  latter  dimension  is 

equal  to  the  number  of  interior  vertices  in  the  triangulation  with  mesh  size  h.  We  then  have 
the  following  relation  between  stiffness  matrices  and  bilinear  forms  and  vector  representations 
of  degrees  of  freedom: 

(Aoflv,w)*a  =  Ao(v,w),  (Acx,r])e2  =  a(x,v)-  (5-8) 

These  identities  hold  for  all  v  G  and  iv  G  V{fR  and  their  vectors  of  degrees  of  freedom, 
v  G  ]Rri'CH  and  w  G  IRncfl,  as  well  as  for  all  x  £  F~onf,  and  ah  V  £  V~onf  and  the  respective 
vectors  of  degrees  of  freedom  x  £  lRnc  and  r]  G  IR™C. 

Since  the  conforming  finite  element  functions  are  in  V£' R .  we  may  expand  each  of  the  canonical 
basis  functions  from  V~oui  via  the  basis  in  V^R.  The  coefficients  in  such  expansions  form 
a  matrix  n  G  mncRXnc  which  represents  the  natural  inclusion  operator  V~oni  C  VffR.  The 
additive  preconditioner  B  :  IRMCii  eA  ]RncH  is  then  defined  in  the  usual  manner: 

B  :=  Dq  1  +  n(Ac’)_1nt.  (5.9) 

Denoting  by  A/o(7^)  the  set  of  interior  vertices  of  7^  we  see  that  the  following  relation  holds  for 
all  p  G  and  p'  G  AfoiT^): 

[A%y  =  a((pp,  pp.)  =  [nTACi?n]p  p,  ,  i.e.  Ac  =  UTACRU.  (5.10) 
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5.1.2.  Definition  of  the  preconditioner:  Operator  notation.  We  first  consider  the  opera¬ 
tor  R corresponding  to  the  Jacobi  smoother.  We  note  that  the  usual  definition  of  the  operator 
associated  with  Jacobi  relaxation  requires  additional  notation,  including  introducing  mass  ma¬ 
trices  and  other  objects,  which  are  not  actually  used  in  the  implementation  of  the  algorithms. 
To  avoid  introducing  these  extra  quantities,  we  will  take  a  slightly  different  path  and  we  define 


R  1AqR  instead  of  R  1 .  We  set: 


R~1Acrv  =  V 


A0(v,fl 


PeR ) 


ee£f 


•Ao((p°R,ipgR) 


„CR  \ 
Ae  = 


[A0Ci?v], 

^  [A^]e„ 


(5.11) 


It  is  clear  from  this  definition  that  the  vector  of  degrees  of  freedom  for  the  function  R~1AqRv 
is  D^A^v,  which  precisely  shows  that  the  matrix  representation  of  R is  Dg  1. 

We  next  consider  the  operator  corresponding  to  the  correction  from  the  conforming  space 
For  its  definition  we  need  to  introduce  the  standard  .Ao(-,  ^-orthogonal  and  1/2(0)- 


prconf 

h 


orthogonal  projections  on  V~°  .  Let  us  denote  these  projections  with  Pc  and  Q  ,  respectively. 
They  are  defined  for  any  v  £  V£R  and  [Pcv\  £  Wconf  and  [Qcv\  £  Wconf  are  as  follows: 


a([Pcv\,rj)  =A0(v,ri),  ([ QCv\,ri /)  =  (v,ij), 
The  operator  form  of  the  preconditioner  is  then 


V  r,eV~ 


conf 


JD  .  T  tCR  _ v  T/C'.R 

B  ■  vh  ^  vh  > 


B  :=  R~l  +  (AcylQc 


(5.12) 


(5.13) 


We  also  have  the  following  identity,  which  follows  directly  from  the  definition  of  Qc  and  Pc : 

QcA$r  =  ACPC.  (5.14) 

Multiplying  from  the  right  the  defining  relation  for  B  in  (5.13)  with  Aq  and  using  the  above 
relation  (5.14)  we  have  that 


BAqH  :=  R~1A%r  +  P 


C 


(5.15) 

From  the  definitions  given  above  it  is  immediate  to  show  that  the  vector  of  degrees  of  freedom 
corresponding  to  the  function  BAq  rv  is  BAq  v,  where  v  is  the  vector  of  degrees  of  freedom 
for  v(x).  This  together  with  (5.8)  shows  that  the  following  identity  holds  for  any  v  £  V^R: 


(BA^v,A^v)4  _  A0(BACrv,v) 


(5.16) 


(Ag^v,  v)^2  A0(v,v) 

Note  that  the  left  side  of  this  identity  uses  the  1 2  inner  product  on  IR"Cfl,  while  calculating  the 
right  side  is  done  through  the  bilinear  form  Aq(-,  •). 


5.1.3.  Spectral  equivalence  results  for  the  two  level  preconditioner.  We  recall  definition 
of  the  set  of  indices  X  of  floating  subdomains  (see  [69])  which  are  the  subdomains  where  the 
coefficient  is  constant  and  which  do  not  touch  the  Dirichlet  boundary: 

X  :=  {  i  :  measd-i(<912  FI  OOfl  =  0  }  .  (5-17) 

The  role  of  the  set  X  will  be  clear  from  the  convergence  analysis  given  later  on  in  this  section. 
The  main  spectral  equivalence  result  for  the  two  level  preconditioner  is  as  follows. 

Theorem  5.2.  Let  mo  =  |X|  be  the  total  number  of  floating  subdomains.  Then,  the  effective 
condition  number  K,eff  (BAq  R)  :=  /Cmo+i(SAg  )  for  the  preconditioner  B  defined  in  (5.13)  is 
uniformly  bounded: 

K.mo+i(BAo )  <  (1  +  |  log  2h/h\)  <C,  forh>h , 
with  C  >  0  independent  of  the  coefficients. 
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The  rest  of  the  section  is  devoted  to  the  proof  of  Theorem  5.2,  which  essentially  involves 
showing  spectral  equivalence  of  the  form 

Aq(v,  v)  <  (B^v,  v )  <  A0(v,  v ). 

It  can  be  shown  (see  Appendix  B,  Lemma  B.2)  that  the  inverse  of  B  satisfies: 

(B~lv,v)  =  inf  [R(v-x,v-x)  +  a(x,x)\,  (5-18) 

x£yconi 

h 

where  R(-,-)  is  the  bilinear  form  associated  with  the  smoother  and  is  defined  as:  lZ(v,w)  := 
(Rv,  w).  In  the  Appendix  B  (see  Lemma  B.l)  it  will  be  shown  that  for  all  v  £  VjfR  we  have 

U(v,  v)=Y^  (¥>?*  <PeR)vi(P  ov,  v)h.  (5.19) 

^£°h 

From  the  above  considerations  it  is  clear  that  the  proof  of  Theorem  5.2  amounts  to  verifying 
the  following  two  conditions. 

(FI)  There  exists  a  constant  c\  such  that  for  all  v  £  r  we  have 

cflA0{v,v)  <  K(v,v). 

(F2)  For  every  v  £  V^R  there  exists  Xv  £  k~conf  (depending  on  v),  such  that 

R(v  ~Xv,v~  Xv)  +  a(xv,  Xv)  <  cqAq{v,  v). 

We  remark  here  that  the  identity  (5.18)  is  a  well  known  result  regarding  two  level  additive 
Schwarz  methods  (cf.  [77,  67,  69]). 

The  validation  of  last  assumption  (F2)  is  more  intricate  and  we  postpone  it  to  the  next 
subsection,  while  the  Lemma  stated  next  shows  that  (FI)  is  true  for  our  choice  of  7 Z(-,  •),  and 
also  gives  an  inequality  which  will  later  be  used  to  verify  (F2). 


Lemma  5.3.  Let  TZ(-,-)  be  the  bilinear  form  defined  via  (5.19).  Then  we  have  the  following 
estimates 


cx  1Aq{v,  v)  <  TZ(v,  v)  and  1 Z(v,v)~h  2|M|oK,  Vu  £  Vfi 


(5.20) 


Proof.  To  show  that  (FI)  holds  we  set  w  =  Yhe  weT>e  and  Cauchy-Schwarz  and  the  arithmetic- 
geometric  inequalities  give 

A0(v,v)  =  ^2^2A0(<PeR,<Pe'R)VeVe> 

e  e' 

<  \/ A^PeR^eR)  \j A0{ip^R,^,R)veVe/ 

e  e' 

<  )  EE  [Al(vfK,  +  A>(k£S,  fe’R)Ve>] 


E  E  MV ?*,  PP)vl  <  E  PeR)vl  =  K(v,  V). 


The  constant  hidden  in  the  above  only  depends  on  the  number  of  neighboring  faces  e! ,  for 
a  given  face  e  £  and  this  constant  is  bounded  by  5  in  2D  and  7  in  3D. 

Notice  that  since  the  mesh  is  quasi-uniform,  for  any  3  v  =  sffiJeveT^R  and  T  £  Th  we 
have 

1MIoa,T-  Ve\\VeR\\o,K,T-  (5-21) 
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On  one  hand,  for  any  basis  function  by  the  inverse  inequality  (see  [29])  we  have  that: 
^2\\t^R\\o  kT  ~  ll^e^llo  K  T-  On  the  other  hand,  directly  calculating  the  norms  gives  us  the 
following  bound  HV^^Hq  K  T  < /i_2||^i?||g  K  T.  Hence, 

^II^IILt-IIV^IILt-  (5.22) 

Therefore  from  the  definition  of  77(-,  •)  given  in  (5.19)  and  the  equivalence  relations  (5.21) 
and  (5.22)  we  get 


Tl(v,  v ) 


=  E  Ve\\^eR\\lK  =  EVe  E  W^VeYo  ,n,T 

e  e  dTDe 

=  E  E  »«2iiw?*ii§,„t  =  e  E 

T&Th  eCdT  TeTh  eC dT 

—  n  ITIIo,k,t  —  n  ITHo.k- 

TeTh 


□ 


5.2.  A  stable  Decomposition.  In  this  subsection  we  show  that  condition  (F2)  is  satisfied 
and  the  main  tool  is  an  operator  Pfc  :  Vjp R  — >•  V~°  ,  that  will  be  shown  to  satisfy  certain 
approximation  and  stability  properties  as  stated  in  the  next  Lemma. 

Lemma  5.4.  There  exists  an  interpolation  operator  Pfc  :  V{fR  — >  V~°  ,  that  satisfies  the 
following  approximation  and  stability  properties: 

Approximation:  ]](/  —  P^)v\\o,K  <  Cah\  log  2/i//t|1,/2||u||i)/lA,  \/vEV/fR,  (5.23) 

Stability:  I-P^Iia  <  Cs\ log 2/j.//?,| 1/2 1 1 u 1 1 VuGl VffR,  (5.24) 

with  constants  Ca  and  Cs  independent  of  the  coefficient  k  and  mesh  size. 

The  definition  of  such  Pjf  and  the  proof  of  the  above  result  is  given  in  the  Appendix  A.  We 
would  like  to  point  out  that  the  operator  Pjf  is  not  used  in  the  actual  implementation  of  the 

preconditioner  B,  as  is  plainly  seen  from  (5.9).  However,  the  operator  Pjf  and  its  approximation 
and  stability  properties  play  crucial  role  in  the  analysis. 

We  now  discuss  and  display  the  main  ingredients  in  the  proof  of  condition  (F2),  with  the 
aid  of  Lemma  (5.4).  Observe  that  on  the  right  hand  side  of  (5.23)  and  (5.24),  the  bounds  are 
given  in  terms  of  the  weighted  full  H 1  -norm  ||f  || i,/j, In  general,  one  cannot  replace  the  norm 
IMIi,h,K  by  the  energy  semi-norm  \v\i^:K,  unless  some  stronger  conditions  on  the  coefficients 
(like  quasi-monotonicity)  are  imposed.  To  be  able  to  replace  the  full  norm  by  the  semi-norm, 
we  introduce  the  subspace  V{fR  C  V^R: 

VfR  :=<«£  V{fR  :  [  vdx  =  0,  Vi  G  X 1  ,  (5.25) 

l  Jtii  J 

where  X  is  the  set  of  indexes  of  the  floating  subdomains  as  defined  in  (5.17).  Recall  that  a  floating 
subdomain  is  a  subdomain  whose  boundary  does  not  intersect  the  Dirichlet  boundary.  The  key 
feature  of  the  subspace  we  have  just  introduced  is  that  the  Poincare- Friedrichs  inequality  for 
the  nonconforming  finite  element  space  [36,  17]  holds  on  each  subdomain,  and  this  allows  us  to 
replace  the  full  norm  ||v|| by  the  energy  semi-norm  |u|i for  any  v  G  VjfR.  We  should 
remark  that  the  condition  on  the  average,  i.e.  fn  vdx  =  0,  is  not  essential;  other  type  of 
conditions  could  be  used  (see  [69])  as  long  as  they  allow  for  the  application  of  a  Poincare- type 
inequality. 


20 


B.  AYUSO  DE  DIOS,  M.  HOLST,  Y.  ZHU,  AND  L.  ZIKATANOV 


At  this  point  we  would  like  to  mention  that  the  dimension  of  VFR  is  related  to  the  number 
of  floating  subdomains  and  in  fact,  dim(V)pR)  =  dim (VffR)  —  mo-  By  restricting  now  the  action 

of  the  operator  P^  to  functions  in  VFR,  we  have  the  following  result,  as  an  easy  corollary  from 
Lemma  5.4. 

Corollary  5.5.  Let  V^R  C  V^R  the  subspace  of  the  Crouziex-Raviart  finite  element  space,  as 
defined  in  (5.25).  Then,  there  exist  an  operator  P p  :  VffR  — >  IL?onf  that  is  stable  in  the  weighted 
H1  -broken  semi-norm: 

|-P/N|ia  <  I  log2h//i|1/2|u|i)/lA  ,  VuG  VfiR  , 
and  satisfies  the  following  approximation  property: 

W  ~  Ph)vb,n<h\\og2h/h\1/2\v\lXK  ,  V  v  eVjfR  . 

Proof.  The  proof  follows  from  Lemma  5.4  together  with  a  standard  application  of  the  Poincare- 
Friedrichs  type  inequalities  inequality  for  nonconforming  functions,  which  allows  for  bounding 
the  norm  ||u||i  ^  K  by  the  semi-norm  |u|i  ^iK.  In  fact,  from  the  definition  of  ||v||i  ^ 

Mi, h,K  =  =  5Z('ellullo,ni  +  'eHi,ft,ni) 

fliCQ  ifl 

+ (^ibiioPi + K\v\\hQf) . 

iex 

For  the  terms  in  both  sums  above,  one  can  apply  Poincare-Friedrichs  type  inequalities.  In  the 
first  sum  above,  for  i  (fX  each  subdomain  touches  the  boundary,  and  hence  || w ||o,f2*  <  Cp\v\i,h,n.r 
For  the  terms  in  the  second  sum,  since  v  G  VFR  are  functions  with  average  zero  on  f we  have 
IMIoa  <  Cp\v\i}h,cii,  i  G  X.  □ 

Remark  5.6.  In  particular,  if  we  consider  the  case  when  the  ratio  h/h  is  a  fixed  constant  (e.g. 
h  =  h )  then  the  approximation  and  stability  properties  are: 

||(I  -  pfov\\0tK  <  Clog  (2)h\v\iihtK  <  h\v\i^K, 

\Phv\i,K  <  Clog(2)|u|i)feA  <  \v\ithjK. 

With  the  aid  of  the  result  in  Corollary  5.5,  we  are  able  to  show  that  the  stability  condition 
(F2)  holds  for  functions  from  the  subspace  VFR. 

Lemma  5.7.  For  any  v  G  VffR  the  stability  condition  (F2)  holds,  namely,  there  exist  a  constant 
Co  >  0  and  a  function  x  =  Ph(v)  £  P~onf  suc h  ^at 

2  f 

K(v  ~  X,v  -  x)  +  o(x,x)  <  c0|u|ph  K,  cq  =  C(1  T  log  -j— )  .  (5.26) 

Proof.  Given  any  v  G  Vp  ,  let  %  G  V~onf  be  defined  as  x  :=  Phv ■  Taking  into  account  the 
approximation  property  of  Pjf  (with  h  =  h )  given  in  (5.5)  together  with  Remark  5.6,  we  have 

R(v  -x,v~x)<  h~2\\v  -  x\\o,k  F  h~2\\v  -  P%v ||qa  <  \v\\Xk  =  A0(v,v). 

The  proof  is  now  complete  by  using  the  stability  of  Pjfv  given  in  Corollary  5.5; 

a(x,x)  =  \Phv\jA  <  \  log2h/h\\v\lAK  =  \log2h/h\A0{v,v). 

□ 

We  have  now  all  ingredients  to  complete  the  proof  of  Theorem  5.2. 


MG  FOR  DG  WITH  JUMP 


21 


Proof  of  Theorem  5.2.  To  estimate  the  maximum  eigenvalue  of  BA qR  let  x  £  V/f  and  v  G  V^R 
be  arbitrary.  We  set  vq  =  (v  —  x),  and  hence  u  =  no  +  x-  Applying  the  Cauchy-Schwarz 
inequality,  the  arithmetic-geometric  inequality,  and  (5.20)  then  yields 

Aq{v,  V )  =  Ao(v0  +  x,  V0  +  x)  =  Ao(^o,  vo)  +  Ao(x,  x)  +  2*4o(n0,  x) 

<  2Ao{vq,vq)  +  2*4.0 (x,  x)  <  ci  (R(v0,  v0)  +  a(x,  x)) , 

Since  x  £  V~onf  was  arbitrary,  taking  the  infimum  with  respect  to  x  and  using  the  definition  of 
(. B~1v,v )  (5.18)  shows  that 

Ao(y,v)  <  (.B_1v,v). 

This  guarantees  that  the  maximum  eigenvalue  Amax(PA(,  )  <  c\  is  uniformly  bounded  inde¬ 
pendently  of  the  coefficient  and  the  mesh  size. 

For  the  lower  bound,  by  Lemma  5.7  we  have  that  for  any  v  G  VfR  with  Xv  =  Phv  we  have 

K(v  -Xv,v-  Xv)  +  |Xvli,«  <  coM l:h,K  =  coAo(v,  v). 

Then  by  the  minimax  principle  (cf.  [52,  §90],  [6,  Lemma  3.13]  or  [46,  Theorem  8.1.2])  and 
noticing  that  dim (V^R)  =  dim(V)pR)  —  mo,  we  have  Amo+i  >  Cq  1.  Therefore,  the  effective 
condition  tCmo+i(BA^R)  can  be  bounded  by 

/CTO0+i  (BA%R)  <  c0ci, 

with  co  given  in  (5.26).  This  completes  the  proof.  □ 

Remark  5.8.  All  the  results  that  we  have  proved  above  also  hold  for  R ~1  defined  by  the  sym¬ 
metric  Gauss-Seidel  method.  This  follows  from  the  fact  that  the  bilinear  forms  77(-,-)  defined 
via  Jacobi  or  symmetric  Gauss-Seidel  method  are  equivalent  with  constants  independent  of  the 
variations  in  the  PDE  coefficient.  For  a  proof  of  such  equivalence  we  refer  to  [71,  Proposi¬ 
tion  6.12]  or  [81,  Lemma  3.3].  Thus,  the  two  level  (and  also  multilevel)  preconditioners  that  use 
Jacobi  or  Gauss-Seidel  as  smoother  are  equivalent  and  hence  all  the  results  here  stated  for  the 
Jacobi  smoother  also  hold  for  the  symmetric  Gauss-Seidel  smoother. 

5.3.  Multilevel  Preconditioner  for  Aq ( • ,  • )  on  VpR  .  We  now  introduce  multilevel  precon¬ 
ditioner,  using  the  two  level  theory  developed  before.  The  multilevel  preconditioner  corresponds 
to  replacing  [Ac]_1  in  (5.13)  with  a  spectrally  equivalent  operator  B(  :  V~oni  *— >  k~onf  corre¬ 
sponding  to  the  additive  BPX  preconditioner  (see  e.g.  [13,  14,  77]). 

The  space  decomposition  that  we  use  to  define  the  multilevel  BPX  preconditioner  is: 

J  J+i 

yCR  =  yCR  +  J2Wj  =  Y,  Wj,  (5.27) 

j=o  j= o 

where  we  have  denoted  Wj  =  V)(conf  (j  =  0, 1,  •  •  •  ,  J)  the  nested  conforming  spaces  with  mesh 
sizes  hj  =  2?h  and  JJj+i  =  V{fR.  Let  us  also  denote  by  A *?  the  operators  corresponding  to  the 
restrictions  of  a(-,  •)  on  Wj  for  j  =  0,  •  •  •  ,  J,  namely 

(. AcfVj,Wj )  =  a(vj,Wj),  Vvj  G  Wj,  Mwj  G  Wj. 

Without  loss  of  generality,  we  assume  that  h  —  2~  J,  and  therefore  |log/i|  ~  J.  Then  the 
operator  form  of  the  multilevel  preconditioner  is: 

J+i 

Bml  ■=  [Aq]  1Qo  +  ^2  Rj  1(5i 

3= 1 


Bml  ■■  V°R  hg  VhCR, 


(5.28) 
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Here  Q3  :  VfiR  ha  Wj  is  the  L2-orthogonal  projection  on  Wj,  j  =  0, , . . ,  J  and  we  set  Qj+i  =  I. 
We  use  exact  solver  on  the  coarsest  grid.  With  all  this  notation  in  hand  one  can  prove  that 


v) 


E 


inf 

j=0  wi  =v 


J+ 1 

a(w0,  w0)  +  E  lZj(wj,Wj ) 

3= 1 


(5.29) 


Here  R-jf,  •),  j  =  1, . . . ,  ( J  +  1)  correspond  to  Jacobi  or  symmetric  Gauss-Seidel  smoother  and 
the  proof  of  (5.29)  is  similar  to  the  proof  in  the  two  level  case. 

Let  us  recall  two  results  that  we  need  later  in  our  convergence  analysis. 


Lemma  5.9  ([78,  Lemma  4.2]).  Let  Rjf,  •)  be  a  Jacobi  or  the  symmetric  Gauss-Siedel  smoother 
for  the  solution  of  the  discretization  (5.6)  on  Wj  space.  Then, 

a(w,  w)  <  Kj(w,w)  <  hj2\\w\\'lK  Wj 

We  also  need  the  following  strengthened  Cauchy  Schwarz  inequality. 


Lemma  5.10  (Strengthened  Cauchy  Schwarz,  cf.  [77,  Lemma  6.2]).  Let  j  =  0, 1,  •  •  •  ,  J  —  1  and 
j  <  l  <  J  there  exists  a  constant  7  G  (0, 1)  such  that 

a(wi,Wj )  <  'yl~j{hf1\\wi\\o,K){h~1\\wj\\0tK),  \/wt  G  Wt,  and  wj  G  Wj.  (5.30) 

The  main  result  of  this  section  is  the  following: 


Theorem  5.11.  Let  Bml  be  the  multilevel  preconditioner  defined  in  (5.28)  and  let  mo  denote 
the  number  of  floating  subdomains  ( the  cardinality  of  the  set  Z  defined  in  (5.17),).  Then,  the 
following  upper  bound  can  be  given  for  the  effective  condition  number  JCmo+i(B mlAqR)  for  the 
multilevel  preconditioner  Bml- 

^mo+l (BMLAqR)  <  C\  loghl2  <  I  log/?]2  , 

where  the  constant  C  >  0  is  independent  of  the  coefficients  and  mesh  size. 

Proof.  To  estimate  the  effective  condition  number,  we  restrict  our  considerations  to  v  G  VPR. 
We  first  give  a  bound  on  the  first  relevant  eigenvalue  \mo+i(Bf^AQR) .  From  (5.29)  we  can  see 
that  to  estimate  this  eigenvalue  we  need  to  find  a  decomposition  v  =  ]Cj=()  w3  w^c^1  stable. 
To  simplify  the  notation  we  set  Pj  :=  P^f  :  V^R  — >  Wj,  for  j  =  0, . . . ,  J,  and  Pj+i  =  I ,  and 
P- 1  =  0.  Then  for  any  v  G  VffR,  we  set 

J+i  J 

v  =  ^(P,  —  Pj-i)v  =  ^2 Wj ,  where  Wj  =  (Pj  —  Pj-\)v. 
j= 0  j= 0 

Clearly,  wj  G  Wj  for  j  =  1,  •  •  •  ,  ( J  +  1)  and  wq  =  Pov  G  Wq.  Now  we  show  that  this  de¬ 
composition  is  stable,  namely  that  exist  Co  >  0  such  that  for  all  v  G  V^R  decomposed  as 
above, 

J  +1 

a(w0,  w0)  +  ^2  R-j{wj-,Wj)  <  CqAq(v,  v)  ,  (5.31) 

3= 1 

To  estimate  TZj+\(wj+\,wj+f)  we  use  (5.20)  from  Lemma  5.3,  together  with  the  approximation 
result  (with  h  =  h)  from  Corollary  5.5  and  Remark  5.6.  Then  by  the  triangle  inequality,  with  the 
approximation  and  stability  estimates  of  Pj  from  Corollary  5.5,  and  the  fact  that  log  hj/h  —  j , 
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gives 


J+ 1 


J+ 1 


a(wo,uio)  +  y~]7 Zj(wj,Wj) 


< 


3= 1 


i=i 

J+i 

<  I^H?A,n  +  2^/r2||u-P.,u| 
i=i 

J+i 


IL 


< 


\lh,K^2\loS^\  <  =  ^o(u,v)  • 

j=i 


Hence,  (5.31)  is  shown  with  Cq  =  J2  ~  |  log/i|2,  and  so  Amo+i(.BMLAQ^)  >  C) 


-2 


We  next  give  a  bound  on  Amax(L>MLAQR).  Let  v  €  W j+ \  =  VGR  be  arbitrary  and  iwj}'t=o  be 

J- 
3= 


any  decomposition  of  v,  namely  v  =  w3  >  with  wj  £  Wj-  By  strengthened  Cauchy-Schwarz 


inequality,  we  have 

(J+ i  J+1  ii  i 

Al0(u,  u)  =  A  -3  ^oK+1>^+1)  +  aKlUo)  +  a 


ILU 


i 


\3= 0  ,?=0 


,j=i  i=i 


j  j 


=  3  A0(wj+i,wj+1)  +  a(w  o,  w0 )  +  EE  a  (uij,  uij) 

\  i=i  j=i 

By  Lemma  5.9,  the  equivalence  (5.20)  and  noticing  that  the  spectral  radius  of  the  matrix 
(7|i_j|)jxj  is  uniformly  bounded  by  (1  —  7)  1  we  get: 

J- 1  J 

A0(v,v)  <  Kj+i(wj+1,wj+i)  +  a(w0,w0)  +  ELtni  (hi  ||uit||o,/t)  {hj  Iloilo, kJ 

i= 0  3=1 

(  3 

<  C1!  K(wj+Uwj+1)  +  a(wo,wo)  +  y^Kj(wj,wj) 

\  3= 1 

Since  the  decomposition  was  arbitrary,  taking  the  infimum  over  all  such  decompositions  together 
with  (5.29)  then  gives 

Aq(v,v)<Ci(B^v,v),  VveVfR, 


which  shows  that  Amax(-BML24-o  )  <  Ci,  and  the  proof  is  complete. 


□ 


Remark  5.12.  Similar  results  hold  also  for  the  multiplicative  multilevel  methods  such  as  the 
U-cycle  and  these  results  can  be  easily  derived  from  estimates  comparing  multiplicative  and 
additive  preconditioners  given  in  [50,  Theorem  4]  or  [27,  Theorem  4.2]. 


6.  Solvers  for  IP-1  Methods 

We  now  introduce  the  different  iterative  methods  for  the  solution  of  (2.5).  As  we  will  see, 
most  of  these  methods  are  based  on  the  methods  constructed  previously  for  the  solution  of  (2.9). 
We  begin  by  describing  the  general  setting  for  the  construction  of  the  solvers.  In  all  cases,  we 
follow  the  ideas  from  [9],  and  we  will  focus  on  the  construction  of  a  preconditioner  (iterator) 
denoted  by  BDG .  For  simplicity  we  consider  the  following  linear  iteration: 

Algorithm  6.1.  Given  initial  guess  uq,  for  k  =  0, 1 . . .  until  convergence: 

1.  Set  ek  =  BDG(f  -  Auk)\ 

2.  Update  uk+ 1  =  uk  +  ek  ■ 
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Here  A  is  the  operator  associated  with  to  the  bilinear  form  of  any  of  the  IP-1  methods. 

We  next  discuss  the  construction  and  properties  of  the  preconditioner  (iterator)  BDG  for 
symmetric  and  non-symmetric  methods. 

6.1.  Solvers  for  the  SIPG  method.  For  the  SIPG  method  (6  =  —  1  in  (2.5)),  BDG  is  used 
in  the  PCG  algorithm  as  a  preconditioner.  From  the  spectral  equivalence  between  A(-,  •)  and 
^4.o ( -,- )  given  in  Lemma  2.2,  it  follows  that  any  of  the  preconditioners  designed  for  Ao(;,  •)  result 
in  an  efficient  solver  for  4(-,-).  In  particular,  due  to  the  block  diagonal  form  of  Ao,  we  focus 
on  block-Jacobi  preconditioners.  We  denote  by  Rz  and  [Rq]  the  operator  corresponding  to  the 
diagonal  of  A{ -,  •)  and  4o(-,  •),  respectively  restricted  to  Zp.  Using  the  decomposition  (3.3)  from 
Proposition  3.1,  we  now  define  the  following  preconditioners: 

Block-Jacobi:  BGC  :=  [R2]”1  +  BQCR  ,  (6.1) 

Block-Jacobi  for  Ao-  BRG  '■=  [-Rq]-1  +  BqQgr  ,.  (6.2) 

Here  QGR  :  VDG  i-A  VGR  is  the  ^-orthogonal  projection  on  VGR.  Here  Bo  =  Bml  as  defined 
in  (5.28)  and  B  refers  to  the  corresponding  multilevel  preconditioner  for  the  symmetric  SIPG-1 
method  (i.e.,  including  the  jump-jump  term). 

The  next  result  is  a  simple  consequence  of  the  analysis  given  in  the  last  two  sections  (Theorems 
5.2  and  5.11)  together  with  Lemma  2.2: 

Theorem  6.2.  Let  BDG  be  the  preconditioner  defined  through  either  (i)  or  (ii).  Let  mo  denote 
the  number  of  floating  subdomains.  Then,  the  following  estimate  holds  for  the  effective  condition 
JCmo+i  (BdgA): 

JCm0+i(BDGA)  <  C\logh\s  <  |log/r|s  , 

where  s  =  0  for  the  two  level  method  and  s  =  2  for  the  midtilevel  method.  The  constant  C  >  0 

above  is  independent  of  the  variation  in  the  coefficients  and  mesh  size. 

6.2.  Solvers  for  the  non-symmetric  IIPG-1  and  NIPG-1  methods.  We  now  discuss  two 
possible  choices  for  the  preconditioner  (iterator)  BDG  for  solving  the  non-symmetric  methods. 
Since  the  operator  notation  is  convenient  in  describing  these  preconditioners,  for  a  given  4(-,  •) 
(corresponding  to  either  IIPG-1  or  NIPG-1  discretization)  we  define  the  operator  A  :  VDG  H > 
VDG  in  a  standard  way: 

( Av,w )  :=  A(v,w),  \/v  G  VDG,  \/w  G  VDG.  (6.3) 

6.2.1.  Preconditioning  with  the  symmetric  part  of  A(-,-).  The  first  preconditioner  we 

consider  is  the  inverse  of  the  symmetric  part  of  A,  and  is  defined  as  follows: 

Bdg  :=  Ag1,  where 

(Asv,w)  :=  ^\A(v,w)  +  A(w,v)],  Vv  G  VDG,  Vw<eVdg. 

We  note  that  from  this  definition  and  (2.15),  we  immediately  have  that  As  is  symmetric  and 
positive  definite  and  hence  defines  a  norm,  which  we  will  refer  to  as  the  Hs'-norm  and  denote  as 
||  •  ||as-  We  briefly  discuss  and  justify  now  why  the  results  given  in  [9]  can  be  extended  to  the 
jump  coefficient  problem  and  henceforth  applied  to  the  present  situation.  We  mainly  focus  on 
the  IIPG  method  since  the  convergence  results  for  NIPG  can  be  deduced  by  little  modification 
of  those  for  IIPG.  We  state  the  following  result  (which  is  an  extension  of  [9,  Theorem  5.1]  to 
the  model  problem  (1.1),  and  show  uniform  convergence  of  the  linear  iteration  in  Algorithm  6.1 
with  iterator  BDG  given  by  (6.4)): 
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Theorem  6.3.  Let  a*  be  a  fixed  value  of  the  penalty  parameter  for  which  the  IIPG-0  bilinear 
form  (2.9)  *4o('o)  is  coercive.  Let  be  the  bilinear  form  of  the  IIPG-1  method  (2.5)  with 

penalty  parameter  a  >  4a*.  Let  BDG  =  be  the  iterator  in  the  linear  iteration  6.1  and  let 
uk  and  uk+\  be  two  consecutive  iterates  obtained  via  this  algorithm.  Then  there  exists  a  positive 
constant  A  <  1  such  that 

||u  -  uk+ i\\DG  <  A||w  -  Uk \\DG  ■  (6.5) 

We  notice  that  the  convergence  is  guaranteed  under  a  technical  (but  mild)  restriction  on  the 
penalty  parameter,  namely  it  should  be  greater  than  4a*. 

A  basic  step  in  the  proof  of  the  above  Theorem  is  to  provide  a  uniform  bound  for  the  skew- 
symmetric  part  of  A(-,  ■),  since  ([BDG]~1u,w)  —  A(u,w)  =  (( As  —  A)u,w)).  This  can  be  done 
proceeding  exactly  as  in  [9,  Lemma  5.6],  only  notational  changes  are  involved.  In  the  end,  the 
result  relies  on  two  ingredients:  a  strengthened  Cauchy-Schwarz  inequality  that  measures  the 
angle  between  the  two  spaces  in  the  decomposition,  namely  Zg  and  VGR,  and  a  simple  Lemma 
that  uses  this  strengthened  Cauchy-Schwarz  inequality  to  give  a  further  bound  for  the  term 
Aq{z,v)  ,  z  €  Zg  ,  v  G  VGR.  This  is  the  step  that  gives  rise  to  the  technical  condition  on  the 
penalty  parameter.  Here  the  key  result  is  the  fact  that  the  functions  in  Zq  are  orthogonal  with 
respect  to  the  inner  product  defined  by  Ao(-,-);  ori  equivalently,  that  the  associated  stiffness 
matrix  hffi  is  diagonal.  This  key  result  is  given  in  Lemma  4.4  for  the  jump-coefficient  problem. 
We  point  out  that  the  proofs  of  the  Cauchy-Schwarz  inequality  and  the  simple  Lemma  mentioned 
above  can  be  carried  out  exactly  as  in  [9,  Lemma  A  &  Lemma  5.1],  respectively,  with  only  small 
changes  in  the  notation.  As  such,  we  omit  these  proofs  here. 

6.2.2.  Block  Jacobi  preconditioner.  The  second  choice  is  simpler  to  implement  and  requires 
less  computational  work  to  apply.  It  is  defined  below,  and  is  naturally  referred  to  as  a  block 
Jacobi  preconditioner. 

Bdg  ■=  Dg1,  where  ... 

(■ Dsv,w )  :=  {Ds(vCR  +  ipz),wCR  +  fiz)  :=  (Agtpz  ,fiz)  +  (AsvCR,wCR). 

In  the  definition  of  Dg  here  we  have  set  v  G  VDG ,  v  =  vGR  +  <pz  and  w  G  VDG ,  w  =  wGR  +  ifz, 
and  vGR,  ipz,  wCR,  and  fiz  are  the  components  from  the  decomposition  (3.3).  As  is  easily  seen 
from  the  definition,  Dg  is  the  block  diagonal  of  Ag.  The  corresponding  matrix  form  of  Dg 
is  denoted  by  Dg.  According  to  the  theoretical  results  stated  and  proved  in  Section  4  and  in 
Section  5,  we  can  solve  a  linear  system  with  respect  to  Dg  efficiently.  We  do  not  present  a 
similar  result  here  for  the  preconditioner  BDG  given  in  (6.6).  Although  plausible,  proving  such 
a  result  does  not  appear  to  be  straightforward;  such  estimates  are  the  subject  of  current  and 
future  research.  We  refer  to  the  numerical  experiments  Section  7  for  further  discussion  and  a 
comparison  of  the  two  preconditioners. 


7.  Numerical  Experiments 


We  consider  the  model  problem  (1.1)  in  the  square  H  =  [—1,  l]2  with  coefficients: 


f  1.0,  Vx  G  [—0.5,  0]2  U  [0, 0.5]2 
(  e,  elsewhere. 


In  all  of  the  following  experiments,  e  varies  from  10-5  up  to  105,  covering  a  wide  range  of  vari¬ 
ations  of  the  coefficients.  In  the  experiments,  we  consider  uniform  refinement  with  a  structured 
initial  triangulation  on  level  0  with  32  elements  and  mesh  size  h  =  2-1.  This  initial  mesh 
resolves  the  jump  interface  of  the  coefficients.  Each  refined  triangulation  is  then  obtained  by 
subdividing  each  element  of  the  previous  level  into  four  congruent  elements.  The  number  of 
degrees  of  freedom  Ng  in  the  DG  discretizations  on  each  level  satisfies  Ng  =  4(Nq  for  t  =  0, 1,  2,  3 
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with  Nq  =  96.  For  convenience  a  list  of  all  tables  that  we  have  used  below  and  the  corresponding 
preconditioners  is  given  in  Table  7.1. 


Table 

Methods 

7.2 

two  level  preconditioner  for  the  CR  problem 

7.3 

Diagonal  preconditioner  for  the  Z  part  in  SIPG  discretization 

7.4 

Diagonal  preconditioner  for  the  Z  part  in  NIPG  discretization 

7.5 

B±u  preconditioner  for  SIPG  discretization 

7.6 

A- norm  of  E  =  (/  —  A)AA)  for  IIPG  with  K  =  1,2 

7.7 

A-norrn  of  E  =  (/  —  AiAA)  for  IIPG  with  K  =  4,  8 

7.8 

A-norm  of  E  =  (I  —  D^1  A)  for  IIPG  with  K  =  4, 8, 16 

7.9 

A-norm  of  E  =  (I  —  Ail1  A)  for  NIPG  with  K  =  1,2 

7.10 

A-norm  of  E  =  (/  —  AjAA)  for  NIPG  with  K  =  4, 8, 16 

7.11 

A-norm  of  E  =  (I  —  D^1  A)  for  NIPG  with  K  =  4,  8, 16 

7.12 

GMRes  with  preconditioner  A^1  for  IIPG 

7.13 

GMRes  with  preconditioner  A^1  for  NIPG 

Table  7.1.  List  of  Tables.  In  the  list  we  refer  to  several  operators:  list  refers 
to  B^°  defined  via  (6.1);  A,  defined  via  (6.3);  A5  and  Dg,  defined  via  (6.4)  and 
(6.6). 


We  have  used  the  basis  (3.4)-(3.5)  for  our  numerical  tests  (for  all  of  the  IP  discretiza¬ 
tions).  Whenever  we  have  used  preconditioned  CG  or  preconditioned  GMRES,  the  iterations 
are  stopped  when  the  initial  residual  is  reduced  by  seven  orders  of  magnitude:  namely,  if  ro  is 
the  initial  residual  and  r ^  is  the  residual  at  iteration  £,  the  iteration  process  (PCG  or  precon¬ 
ditioned  GMRES)  is  terminated  at  iteration  k  if  ||  IU2  / 1 1 7'°  1 1  ^2  <  10-7.  For  the  non-symmetric 
discrete  schemes  (NIPG  and  IIPG),  we  have  used  GMRES  with  restart  at  every  10  iterations, 
and  the  maximum  number  of  iterations  is  set  to  30.  The  experiments  were  carried  out  on  an 
IMAC  (OS  X)  with  2.93  GHz  Intel  Core  i7,  and  8  GB  1333  MHz  DDR3. 

7.1.  Solver  for  IP(/3)-0  method.  We  first  consider  IP (/3)-0  method.  For  this  set  of  experi¬ 
ments  we  have  set  the  penalty  parameter  a  =  8.  We  use  algorithm  4.1  given  in  Section  4  to  solve 
the  linear  system  arising  from  the  IP(/5)-0  discretization.  Due  to  the  block  structure(4.1)  of  Ao 
(matrix  representation  of  Ao  in  the  basis  (3.4)-(3.5))  we  only  need  to  numerically  verify  the 
effectiveness  of  the  solvers  for  each  block;  Aqv  and  Aq2.  Recall  that  for  any  choice  of  8  =  0,  ±1, 
the  block  A'^v  is  the  same  (since  it  is  the  stiffness  matrix  of  the  Crouzeix-Raviart  discretization 
(5.1)),  while  the  block  Aq2  is  different  for  different  values  of  9. 

The  system  Aq11  arising  from  the  restriction  of  Ao(-,  •)  to  the  Crouziex-Raviart  space  is  solved 
by  a  PCG  algorithm  with  the  two  level  preconditioner  defined  in  (5.9).  For  the  two  level 
preconditioner  we  use  two  symmetric  Gauss-Seidel  steps  as  smoother.  In  Table  7.2  we  re¬ 
port  the  estimated  condition  number  /C(BAqv)  and  the  effective  condition  number  (denoted  by 
/Ci(BAqw)).  Observe  that  the  estimated  condition  number  /C(BAq^)  deteriorates  with  respect 
to  the  magnitude  of  the  jump  in  coefficient.  In  contrast,  and  as  predicted  by  our  theory,  the 
effective  condition  number  /^(BAq1’)  is  uniformly  bounded  with  respect  to  both  the  mesh  size 
and  the  jump  of  the  coefficient,  as  predicted  by  Theorem  5.2.  To  explain  the  deterioration  of 
/C(BAqi;)  with  respect  to  the  jump  in  the  coefficient,  we  have  shown  in  Figure  7.1,  the  spectrum 
of  the  preconditioned  system  for  e  =  ICE5  and  the  mesh  size  h  =  2-5.  Note  that  there  is  only 
one  (very  small)  eigenvalue  close  to  zero  (which  may  be  related  to  the  fact  that  there  are  only  2 
different  values  for  the  coefficients).  The  systems  corresponding  to  Aq2  are  solved  by  a  PCG  al- 
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Figure  7.1.  Eigenvalue  distribution  of  BA™  for  e  =  1CT5  and  h  =  2-5. 

gorithm  using  its  diagonal  as  a  preconditioner.  The  estimated  condition  numbers  of  B71Aq2 
for  Type-0  SIPG(/3)  and  NIPG(/3)  discretizations  are  reported  in  Table  7.3  and  Table  7.4,  re¬ 
spectively.  Observe  that  the  condition  numbers  of  DJ1Aq2  are  uniformly  bounded  and  close  to 
1,  which  confirms  the  result  established  in  Lemma  4.3;  i.e.,  that  Aq2  is  spectrally  equivalent  to 
its  diagonal. 

7.2.  Solvers  for  SIPG(/3)-l.  We  now  study  the  effectiveness  and  robustness  of  the  proposed 
solvers  for  the  SIPG-1  discretization.  As  before,  we  have  set  the  penalty  parameter  ae  =  8  in 
(2.5)  for  this  set  of  experiments.  We  consider  the  B\)(’  defined  via  (6.1)  in  Section  6.1.  The 
matrix  form  of  this  preconditioner  is  denoted  Bf^. 

Table  7.5  gives  the  estimated  condition  numbers  of  /C(BPGA)  and  the  number  of  PCG  iter¬ 
ations  required  for  convergence.  As  can  be  seen  from  these  two  tables,  the  condition  numbers 
of  the  preconditioned  system  deteriorate  rapidly  when  e  becomes  smaller.  In  the  same  table 
we  give  the  effective  condition  numbers  /Ci(Bf>tTA),  and  observe  that  the  effective  condition 
numbers  are  nearly  uniformly  bounded  with  respect  to  the  coefficients  and  mesh  size.  These 
results  verify  the  theory  predicted  by  Theorem  6.2. 

7.3.  Solvers  for  Nonsymmetric  IP(/3)-l  Methods.  We  now  present  some  numerical  tests 
for  the  nonsymmetric  IP(/3)-l  methods:  IIPG(/3)-l  and  NIPG(/3)-l.  We  consider  the  linear 
iteration  given  in  6.1  with  the  iterator  BDG  as  defined  in  (6.4).  As  we  shall  see,  the  perfor¬ 
mance  of  the  preconditioners  depends  on  the  value  of  the  penalty  parameter.  All  the  tests  are 
computed  with  a  =  I\a*,  with  a*  given  at  the  beginning  of  §7.3.  The  value  of  a*  is  close  to 
the  smallest  value  required  for  ensuring  positive  semi-definiteness  of  the  symmetric  part  of  A. 
In  the  experiments,  we  take  K  =  1,  2, 4, 8, 16. 

Since  our  examples  involve  a  fixed  domain  (unit  square),  the  value  of  a*  can  be  well  approxi¬ 
mated.  For  the  IIPG(/3)-l  discretization  with  mesh  size  h  =  2_1  we  take  a*  =  0.9.  For  all  other 
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(a)  e  <  1 


e 

levels 

0 

1 

2 

3 

h 

2-i 

2"2 

2"3 

2-4 

10"5 

A  (BA™) 
/Ci  (BA™) 

5.49e+04  (18) 

4.71 

4.65e+04  (24) 
4.01 

3.79e+04  (26) 
3.6 

3.24e+04  (26) 
3.47 

10"4 

A  (BA™) 
Ai(BA™) 

5.49e+03  (16) 
4.71 

4.66e+03  (23) 
4.01 

3.79e+03  (23) 
3.6 

3.25e+03  (24) 
3.47 

10"3 

A  (BAS”) 
Ai(BA™) 

551  (15) 

4.69 

469  (20) 

4 

383  (21) 

3.6 

328  (21) 

3.47 

10"2 

A  (BA™) 
Ai(BA™) 

57.2  (14) 

4.52 

49.5  (18) 

3.95 

41.2  (19) 
3.57 

36.2  (19) 

3.44 

10-1 

A  (BA™) 
Ai(BA™) 

7.63  (13) 

3.48 

7.34  (16) 
3.54 

6.71  (17) 

3.32 

6.47  (17) 

3.27 

(b)  e  >  1 


e 

levels 

0 

1 

2 

3 

h 

2_i 

2"2 

2"3 

2-4 

i 

A  (BA™) 
Ai(BA™) 

2.58  (11) 
2.36 

2.77  (12) 
2.56 

2.79  (13) 
2.56 

2.94  (14) 
2.69 

101 

A  (BA™) 
Ai(BA™) 

3.13  (10) 
2.16 

4.62  (15) 
3.51 

4.85  (16) 
3.41 

5.06  (16) 
3.4 

102 

A  (BA™) 
Ai(BA™) 

3.73  (11) 
3.36 

6.06  (16) 
3.8 

7.13  (17) 
3.68 

8.3  (18) 
3.66 

103 

E{BA^V) 

Ai(BAgU) 

3.8  (11) 
3.39 

6.3  (16) 
3.84 

7.55  (18) 
3.71 

9.02  (18) 
3.69 

104 

A(BAgu) 

Ai(BA™) 

3.8  (11) 
3.4 

6.32  (16) 
3.84 

7.59  (18) 
3.71 

9.09  (19) 
3.69 

105 

A(BAgu) 

Ai(BAgU) 

3.81  (11) 
3.4 

6.33  (16) 
3.85 

7.6  (18) 
3.71 

9.1  (20) 
3.69 

Table  7.2.  Estimated  condition  numbers  A(BA™)  (number  of  PCG  iterations) 
and  effective  condition  numbers  Ai(BA™)  for  the  block  A™  in  algorithm  4.1. 


mesh  sizes  that  we  use  in  our  numerical  tests  below  (e.g.  h  =  2  2,  h  =  2  3  or  h  =  2  4)  we  take 
a*  =  1.3. 

To  verify  the  theory  for  Algorithm  6.1  with  the  symmetric  part  of  A  as  an  iterator,  (6.4),  we 
have  computed  the  A-norrn  of  the  error  propagation  operator.  We  denote  this  operator  here 
with  E  and  define  it  in  a  usual  manner: 

E  =  I  —  BdgA  =  I  -  A],1  A. 

Below,  we  have  tabulated  ||.E||^  for  several  values  of  the  parameters  of  interest  (e.g.  e,  mesh 
size,  penalty  parameter).  This  norm  gives  us  the  contraction  number  of  the  linear  iteration  in 
Algorithm  6.1.  That  is,  the  estimate  in  Theorem  6.3  holds  with  A  =  ||.E||^  .  This  norm  is 
computed  as  the  maximum  eigenvalue  of  the  generalized  eigenvalue  problem  given  below: 

(I  -  BdgA)tAs(I  -  BdgA)u  =  A Asu, 

where  the  corresponding  definition  of  the  operators  A ,  As,  Ds  and  the  preconditioners  BDG 
was  given  earlier  in  (6.3),  (6.4),  and  (6.6)  respectively.  For  the  IIPG-1  method,  the  Ag-norms 
of  the  error  propagation  operator  with  preconditioner  (6.4)  and  (6.6)  are  given  in  Tables  7.6 
and  in  Table  7.7,  for  the  different  values  of  K.  Notice  that  for  K  =  1  (Table  7.6)  the  iteration 
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(a)  e  <  1 


levels 

h 

€ 

ict3 

10-4 

10-3 

10~2 

lO"1 

0 

2_i 

1.73  (14) 

1.73  (14) 

1.73  (12) 

1.73  (11) 

1.73  (10) 

1 

2-2 

1.72  (15) 

1.72  (14) 

1.72  (13) 

1.72  (12) 

1.72  (10) 

2 

2-3 

1.72  (15) 

1.72  (14) 

1.72  (13) 

1.72  (11) 

1.72  (10) 

3 

2-4 

1.72  (15) 

1.72  (13) 

1.72  (12) 

1.72  (11) 

1.71  (10) 

(b)  e  >  1 


levels 

h 

e 

i 

101 

102 

103 

104 

10° 

0 

2_1 

1.73  (9) 

1.72  (10) 

1.73  (11) 

1.73  (12) 

1.73  (13) 

1.73  (13) 

1 

2-2 

1.72  (10) 

1.72  (10) 

1.72  (11) 

1.72  (12) 

1.72  (13) 

1.72  (14) 

2 

2"3 

1.71  (10) 

1.7  (10) 

1.72  (11) 

1.71  (12) 

1.72  (14) 

1.72  (15) 

3 

2-4 

1.71  (10) 

1.69  (10) 

1.69  (11) 

1.69  (12) 

1.7  (14) 

1.69  (16) 

Table  7.3.  Estimated  condition  numbers  /C(DZ  1Aqz)  (number  of  PCG  itera¬ 
tions)  for  the  block  Aq2  in  Type-0  SIPG(/3)  discretization. 


(a)  e  <  1 


levels 

h 

e 

10"5 

10"4 

10"3 

10"2 

lO”1 

0 

2_1 

1.4  (12) 

1.4  (12) 

1.4  (10) 

1.4  (10) 

1.4  (8) 

1 

2~2 

1.4  (12) 

1.4  (11) 

1.4  (10) 

1.4  (10) 

1.4  (9) 

2 

2"3 

1.4  (12) 

1.4  (11) 

1.4  (10) 

1.39  (9) 

1.39  (9) 

3 

2-4 

1.4  (12) 

1.39  (11) 

1.39  (10) 

1.39  (9) 

1.38  (8) 

(b)  e  >  1 


levels 

h 

e 

i 

101 

102 

103 

104 

10s 

0 

2_1 

1.4  (7) 

1.4  (8) 

1.4  (9) 

1.4  (10) 

1.4  (11) 

1.4  (11) 

1 

2"2 

1.39  (8) 

1.38  (8) 

1.38  (9) 

1.38  (10) 

1.4  (11) 

1.4  (12) 

2 

2"3 

1.39  (8) 

1.37  (8) 

1.36  (9) 

1.35  (10) 

1.35  (12) 

1.35  (13) 

3 

2-4 

1.38  (8) 

1.36  (8) 

1.35  (9) 

1.34  (10) 

1.35  (12) 

1.35  (13) 

Table  7.4.  Estimated  condition  numbers  /C( D,  xAq2)  (number  of  PCG  itera¬ 
tions)  for  the  block  Aq2  in  Type-0  NIPG(/3)  discretization. 


with  BDa  =  Ag1  is  not  convergent  in  Ag-norm.  For  I\  >  4  (as  required  by  our  theory  for  As) 
we  have  a  convergent  algorithm.  As  expected,  (6.4)  converges  faster  for  larger  a. 

The  same  set  of  experiments  is  performed  for  the  NIPG  discretization.  The  estimates  for 
the  A-norrn  of  the  error  propagation  operator  corresponding  to  the  iterator  (6.4)  are  given  in 
Tables  7.9  and  7.10.  As  can  be  observed  from  these  numerical  results,  a  larger  value  of  K  (and  so 
of  the  penalty  parameter)  is  needed  than  the  one  for  the  IIPG  in  order  to  produce  a  convergent 
linear  iteration.  From  the  results  reported  in  the  tables,  it  can  be  observed  that  (provided 
K  >  4  for  IIPG  and  K  >  16  for  NIPG)  the  number  of  iterations  required  for  convergence  does 
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(a)  e  <  1 


6 

levels 

0 

1 

2 

3 

h 

2-i 

2-2 

2-3 

2-4 

lcr5 

JC( BfGA) 
Ai(BfGA) 

2.85e+04  (44) 
6.27 

3.37e+04  (44) 
6.33 

3.1e+04  (46) 
6.45 

2.85e+04  (46) 
6.49 

io-4 

A(BGGA) 

Ai(BfGA) 

2.85e+03  (37) 
6.26 

3.37e+03  (38) 
6.32 

3.1e+03  (39) 
6.45 

2.86e+03  (40) 
6.48 

icr3 

A(BfGA) 

/Ci(BfGA) 

288  (33) 

6.24 

340  (34) 

6.3 

313  (34) 
6.42 

289  (32) 

6.46 

lcr2 

A(BfGA) 

Ai(BfGA) 

32  (27) 

6.08 

37  (28) 

6.1 

34.5  (28) 
6.21 

32.3  (27) 

6.25 

KT1 

£(BGGA) 

Ai(BfGA) 

7.25  (22) 

5.62 

7.33  (22) 

5.6 

7.21  (22) 
5.71 

7.13  (22) 

5.73 

(b)  e  >  1 


e 

levels 

0 

1 

2 

3 

h 

2_i 

2~'z 

2~6 

2-4 

i 

£(BGGA) 

Ai(BfGA) 

5.53  (19) 
5.17 

5.76  (20) 
5.45 

5.8  (20) 
5.46 

5.83  (20) 
5.46 

101 

A(BfGA) 

Ai(BfGA) 

6.66  (22) 
5.91 

7.16  (23) 
6.2 

7.16  (23) 
6.25 

7.43  (23) 
6.27 

102 

A(BGGA) 

£i(BfGA) 

6.88  (25) 
6.32 

8.65  (27) 
6.48 

10.3  (27) 
6.55 

12  (28) 
6.55 

103 

A(BfGA) 

Ai(BfGA) 

6.38  (27) 
5.51 

8.98  (30) 
6.53 

11.1  (31) 
6.59 

13.5  (32) 
6.59 

104 

A(BfGA) 

Ai(BfGA) 

6.91  (30) 
6.38 

9.02  (33) 
6.54 

11.2  (35) 
6.6 

13.7  (36) 
6.59 

105 

£(BGGA) 

Ai(BfGA) 

6.91  (33) 
6.38 

9.02  (36) 
6.54 

11.3  (39) 
6.6 

13.8  (40) 
6.59 

Table  7.5.  Estimated  condition  number  /C(B^  A)  (number  of  PCG  iterations) 
and  the  effective  condition  number  /Ci(BPgA). 


not  grow  with  respect  to  the  mesh  size  or  the  coefficient  jump  and  that  E  =  I  —  A~^lA  is  uniform 
contraction  for  such  values  of  K. 

Similar  results  hold  for  the  block  Jacobi  preconditioner  BDG  =  D^1  described  in§6.2.2.  From 
the  rates  presented  in  Table  7.8  and  in  Table  7.11  one  may  conclude  that  for  the  linear  iteration 
with  the  block  Jacobi  preconditioner  D$,  the  error  propagation  operator  E  =  I  —  D^1  A  is  a 
uniform  contraction  in  the  Ag-norm.  The  numerical  tests  indicate  that  such  estimates  on  the 
rate  of  convergence  hold  for  sufficiently  large  values  of  K  (same  values  as  for  the  preconditioning 
with  the  symmetric  part  As)  and  for  both  IIPG  and  NIPG  discretizations. 

Once  we  have  numerically  verified  that  the  linear  iteration  converges  (provided  K  >  4  for 
IIPG  and  K  >  16  for  NIPG),  we  test  the  use  of  A^1  as  a  GMRES  preconditioner.  We  have 
shown  the  number  of  GMRES  iterations  required  for  convergence  for  different  values  of  I\  with 
the  preconditioner  given  by  (6.4)  in  table  7.12.  The  symbol  x  in  Table  7.12  means  that  GMRES 
fail  to  converge  with  these  parameters  for  >  30  iterations.  It  is  clearly  seen  that  preconditioned 
GMRES  converges  uniformly  with  respect  to  the  mesh  size  and  robust  with  respect  to  the  jump 
in  the  coefficient,  provided  that  K  >  4  for  IIPG  and  K  >  16  for  NIPG. 
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(a)  ||£||is  for  IIPG:  ae  =  a* 


levels 

h 

e 

10-° 

10-4 

10"3 

10-2 

10”1 

1 

101 

102 

103 

104 

10b 

0 

2_1 

2.5 

2.5 

2.5 

2.5 

2.5 

2.5 

2.4 

2.4 

2.4 

2.4 

2.4 

1 

2-2 

1.0 

1.0 

1.0 

1.0 

1.0 

1.1 

1.0 

1.0 

1.0 

1.0 

1.0 

2 

2"3 

1.2 

1.2 

1.2 

1.2 

1.2 

1.2 

1.2 

1.2 

1.2 

1.2 

1.2 

(b)  ||£|fts  for  IIPG:  ae  =  2a* 


levels 

h 

e 

10"5 

lO"4 

10"3 

10~2 

10"1 

1 

101 

102 

103 

104 

105 

0 

2-i 

0.53 

0.53 

0.53 

0.53 

0.53 

0.53 

0.51 

0.52 

0.52 

0.52 

0.52 

1 

2"2 

0.33 

0.33 

0.33 

0.33 

0.34 

0.34 

0.33 

0.34 

0.34 

0.34 

0.34 

2 

2"3 

0.37 

0.37 

0.37 

0.37 

0.37 

0.37 

0.36 

0.37 

0.37 

0.37 

0.37 

3 

2-4 

0.39 

0.39 

0.39 

0.39 

0.39 

0.39 

0.38 

0.39 

0.39 

0.39 

0.39 

Table  7.6.  Norm  of  the  error  propagator  E  =  (I  —  for  A  corresponding 

to  IIPG  discretization,  with:  (a)  ae  =  a*  and  (b)  ae  =  2a*. 


(a)  ||£||is  for  IIPG:  ae  =  4a* 


levels 

h 

e 

10"5 

lO"4 

10"3 

io-2 

IO"1 

1 

101 

102 

103 

104 

105 

0 

2_1 

0.20 

0.20 

0.20 

0.20 

0.20 

0.20 

0.19 

0.19 

0.19 

0.19 

0.19 

1 

2"2 

0.14 

0.14 

0.14 

0.14 

0.14 

0.14 

0.14 

0.14 

0.14 

0.14 

0.14 

2 

2"3 

0.16 

0.16 

0.16 

0.16 

0.16 

0.15 

0.15 

0.16 

0.16 

0.16 

0.16 

3 

2-4 

0.16 

0.16 

0.16 

0.16 

0.16 

0.16 

0.16 

0.16 

0.16 

0.16 

0.16 

(b)  \\E\\2As  for  IIPG:  ae  =  8a* 


levels 

h 

e 

10"5 

IO"4 

10"3 

io-2 

io-1 

1 

101 

102 

103 

104 

105 

0 

2_i 

0.09 

0.09 

0.09 

0.09 

0.09 

0.09 

0.09 

0.09 

0.09 

0.09 

0.09 

1 

2"2 

0.07 

0.07 

0.07 

0.07 

0.07 

0.07 

0.07 

0.07 

0.07 

0.07 

0.07 

2 

2"3 

0.07 

0.07 

0.07 

0.07 

0.07 

0.07 

0.07 

0.07 

0.07 

0.07 

0.07 

3 

2-4 

0.08 

0.08 

0.08 

0.08 

0.08 

0.07 

0.07 

0.07 

0.07 

0.07 

0.07 

Table  7.7.  Norm  of  the  error  propagator  E  =  (I  —  A^A)  for  A  corresponding 
to  IIPG  discretization,  with:  (a)  ae  =  4a*  and  (b)  ae  =  8a*. 


We  conduct  the  same  set  of  experiments  for  NIPG(/3)  discretization.  The  number  of  GMRES 
iterations  required  for  convergence  for  different  values  of  K  with  the  preconditioner  (6.4)  is 
shown  in  Table  7.13.  Clearly  with  such  a  preconditioner  the  GMRES  method  is  uniformly 
convergent  with  respect  to  the  problem  parameters  for  K  >  16. 
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(a)  \\E\\2As  for  IIPG:  ae  =  4a* 


levels 

h 

e 

10-s 

10"4 

10"3 

10~'z 

lO"1 

1 

101 

IT 

103 

104 

10° 

0 

2-i 

0.81 

0.81 

0.81 

0.81 

0.80 

0.80 

0.80 

0.80 

0.80 

0.80 

0.80 

1 

2~z 

0.67 

0.67 

0.67 

0.67 

0.67 

0.67 

0.67 

0.67 

0.67 

0.67 

0.67 

2 

2-3 

0.68 

0.68 

0.68 

0.68 

0.68 

0.68 

0.68 

0.68 

0.68 

0.68 

0.68 

3 

2-4 

0.68 

0.68 

0.68 

0.68 

0.68 

0.68 

0.68 

0.68 

0.68 

0.68 

0.68 

(b)  \\E\\2As  for  IIPG:  ae  =  8a* 


levels 

h 

e 

10"s 

10-4 

10"3 

io-'z 

10"1 

1 

101 

IT 

103 

10* 

10° 

0 

2-1 

0.59 

0.59 

0.59 

0.59 

0.59 

0.58 

0.58 

0.58 

0.59 

0.59 

0.59 

1 

2~z 

0.55 

0.55 

0.55 

0.55 

0.55 

0.55 

0.55 

0.55 

0.55 

0.55 

0.55 

2 

2~3 

0.56 

0.56 

0.56 

0.56 

0.56 

0.56 

0.56 

0.56 

0.56 

0.56 

0.56 

3 

2-4 

0.56 

0.56 

0.56 

0.56 

0.56 

0.56 

0.56 

0.56 

0.56 

0.56 

0.56 

(c)  \\E\\2As  for  IIPG:  ae  =  16a* 


levels 

h 

e 

10“° 

10-4 

10"3 

10~'z 

io-1 

1 

101 

10z 

103 

104 

10° 

0 

2-1 

0.52 

0.52 

0.52 

0.52 

0.51 

0.50 

0.50 

0.50 

0.50 

0.50 

0.50 

1 

2~z 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

0.50 

0.50 

0.50 

0.50 

0.50 

2 

2"3 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

3 

2-4 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

0.51 

Table  7.8.  Norm  of  the  error  propagator  E  =  (/  —  DS1A)  for  A  corresponding 
to  IIPG  discretization,  with:  (a)  ae  =  4a*,  (b)  ae  =  8a*  and  (c)  ae  =  16 a*.  . 


Appendix  A.  Construction  of  the  transfer  operator  Pft. 

We  now  give  the  definition  of  an  operator  P[‘  that  will  be  shown  to  satisfy  the  required 
approximation  and  stability  properties  (5.23)-(5.24).  We  start  by  introducing  some  notation 
required  for  our  construction.  For  a  given  conforming  triangulation  Th,  we  denote  by  N(Th)  the 
set  of  vertices  of  the  partition  Th,  and  by  C(Th)  the  set  of  barycenters  of  the  elements  T  £  Th- 
We  still  denote  by  £}x  the  set  of  edges/faces  of  Th-  For  each  vertex  p  £  A f,  let  up  :=  |^J  T  and 

T3p 

lot  =  [J  ojp  for  each  T  £  Th-  Similarly,  on  the  interface  Op  and  Oe  denote,  respectively,  the 

pST 

local  patches  associated  with  the  vertex  p  £  M  and  the  edge/face  e  £  £h  on  the  interface.  We 
now  start  building  the  operator  Pfr-  it  is  constructed  in  several  steps  as  the  composition  of  a 
particular  inclusion  operator  and  the  Scott.-Zhang  quasi-interpolation  operator.  The  first  basic 
idea  is  to  embed  VPR  into  a  higher  order  conforming  finite  element  space  on  the  same  mesh 
Th-  Following  [17],  we  consider  the  space  of  piecewise  quadratic  polynomials  on  Th,  which  we 
denote  by  V^ouf’2.  To  be  able  to  use  the  results  in  [17]  for  the  jump-coefficients  problem,  we 
consider  this  inclusion  at  the  subdomain  level.  Let  Et  :  V^R{£li)  >  V^on  be  the  inclusion 
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(a)  ||£||is  for  NIPG:  ae  =  a* 


levels 

h 

€ 

10“° 

to-'4 

10"J 

10~2 

10-1 

1 

101 

10J 

10^ 

104 

10b 

0 

2~2 

2.7 

2.7 

2.7 

2.7 

2.7 

2.7 

2.7 

2.7 

2.7 

2.7 

2.7 

1 

2~d 

2.0 

2.0 

2.0 

2.0 

2.0 

2.0 

2.0 

2.0 

2.0 

2.0 

2.0 

2 

2-4 

2.2 

2.1 

2.1 

2.2 

2.2 

2.1 

2.1 

2.1 

2.2 

2.2 

2.2 

(b)  \\E\\2As  for  NIPG:  ae  =  2a * 


levels 

h 

e 

10"5 

IO-4 

10-a 

10~2 

10-1 

1 

101 

102 

10a 

104 

105 

0 

2_i 

1.3 

1.3 

1.3 

1.3 

1.3 

1.3 

1.3 

1.3 

1.3 

1.3 

1.3 

1 

2~2 

0.98 

0.98 

0.98 

0.98 

0.99 

0.99 

0.97 

0.98 

0.99 

0.99 

0.99 

2 

2_a 

1.1 

1.1 

1.1 

1.1 

1.1 

1.1 

1.1 

1.1 

1.1 

1.1 

1.1 

3 

2-4 

1.1 

1.1 

1.1 

1.1 

1.1 

1.1 

1.1 

1.1 

1.1 

1.1 

1.1 

Table  7.9.  Norm  of  the  error  propagator  E  =  (I  —  A^A)  for  A  corresponding 
to  NIPG  discretization,  with:  (a)  ote  =  a*  and  (b)  ae  =  2ot* . 


operator  defined  on  each  subdomain  by: 

Ei[v(p)]  =  V  vr(p)  ypeAf  ( Th(Sli ))  u  C  (Th(fii))  •  (A.l) 

The  above  definition  applies  for  all  interior  points  p  £  M  (7/?. (0^))  U  C  (7 h(^i))-  We  introduce 
the  set  E‘p  :=  {T  £  7/,  :  dT  and  meas^-i  (dT  n  cftlj)  >  0}  and  we  note  that  the  cross-points 

do  not  belong  to  this  set.  At  the  boundary  vertices  p  £  dZll ,  we  define  Ei  in  the  following  way: 

(Eiv)(p)  =  ^  vT{p)  Vp  £  AA  (Th{dfli))  U  C  (Th(dfli)) . 

Observe  that,  with  this  construction,  we  have  that  rji  =  EiV  £  V^onf’2(Hj)  is  conforming  at  each 
subdomain.  However,  the  global  function  77^  :=  Pi  will  generally  be  multi-valued  at  the  cross 
points.  The  modification  of  the  definition  of  Et  at  boundary  points  is  done  to  guarantee  that 
the  global  function  7j  will  be  at  least  conforming  every  non-cross  point  of  the  partition  and  so, 
in  particular,  in  the  interior  of  the  boundary  of  each  subdomain.  The  next  result  is  taken  from 
[17],  where  the  author  studies  uniform  preconditioners  and  solvers  for  the  CR  discretization  of 
the  Poisson  problem.  It  guarantees  the  stability  and  approximation  properties  of  the  inclusion 
operator  Ei.  The  proof  can  be  found  in  [17]: 

Lemma  A.l.  There  exists  a  linear  operator  Ei  :  V^R(Tli)  :— >  V^onf’2(Hj)  defined  on  each 
subdomain  such  that  for  any  v  £  V^R 

(i)  I Eiv\1>n.  ~  Mi ,h,su 

(ii)  \\v-EiV\ |0A  <h  \v\lxn. 

We  note  that  VTonf  c  V.con  ’2  and  therefore  a  standard  restriction  operator  would  allow  us 
h  n 

to  represent  the  constructed  p  as  a  function  in  V~oni  (except  at  cross-points).  However,  to 
guarantee  that  the  constructed  operator  satisfies  (5.23)-(5.24),  we  use  the  Scott-Zhang  quasi¬ 
interpolation  operator ,  again  at  the  subdomain  level  and  also  at  the  interior  of  the  interfaces,  say 
T  =  dQi  n  dTlj  i  j,  between  any  pair  of  subdomains.  We  denote  by  Qi  :  — »  V~onf(Hi) 
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(a)  \\E\\2As  for  NIPG:  ae  =  4a* 


levels 

h 

e 

10"3 

10"4 

10"3 

10“z 

lO"1 

1 

101 

103 

103 

104 

10° 

0 

2-i 

0.66 

0.66 

0.66 

0.66 

0.66 

0.66 

0.64 

0.64 

0.64 

0.64 

0.64 

1 

2~z 

0.49 

0.49 

0.49 

0.49 

0.49 

0.49 

0.48 

0.49 

0.49 

0.49 

0.49 

2 

2-3 

0.54 

0.54 

0.54 

0.54 

0.54 

0.53 

0.52 

0.53 

0.54 

0.54 

0.54 

3 

2-4 

0.56 

0.56 

0.56 

0.56 

0.56 

0.55 

0.55 

0.56 

0.56 

0.56 

0.56 

(b)  \\E\\as  for  NIPG:  ae  =  8a* 


levels 

h 

e 

10"3 

10-4 

10"3 

io-'z 

10"1 

1 

101 

10z 

103 

104 

103 

0 

2_1 

0.32 

0.32 

0.32 

0.32 

0.32 

0.32 

0.31 

0.31 

0.31 

0.31 

0.31 

1 

2~z 

0.24 

0.24 

0.24 

0.24 

0.24 

0.24 

0.24 

0.24 

0.24 

0.24 

0.24 

2 

2-3 

0.27 

0.27 

0.27 

0.27 

0.27 

0.26 

0.26 

0.27 

0.27 

0.27 

0.27 

3 

2-4 

0.28 

0.28 

0.28 

0.28 

0.28 

0.28 

0.27 

0.28 

0.28 

0.28 

0.28 

(c)  ||£||is  for  NIPG:  ae  =  2a* 


levels 

h 

e 

10”° 

10-4 

10"3 

10-z 

to-1 

1 

101 

10z 

103 

104 

10° 

0 

2_1 

0.16 

0.16 

0.16 

0.16 

0.16 

0.16 

0.15 

0.15 

0.15 

0.15 

0.15 

1 

2~z 

0.12 

0.12 

0.12 

0.12 

0.12 

0.12 

0.12 

0.12 

0.12 

0.12 

0.12 

2 

2"3 

0.14 

0.14 

0.14 

0.14 

0.14 

0.13 

0.13 

0.13 

0.13 

0.13 

0.13 

3 

2-4 

0.14 

0.14 

0.14 

0.14 

0.14 

0.14 

0.14 

0.14 

0.14 

0.14 

0.14 

Table  7.10.  Norm  of  the  error  propagator  E  =  (I  —  Asl  A)  for  A  corresponding 
to  NIPG  discretization,  with:  (a)  ae  =  4a*;  (b)  ae  =  8a*;  and  (c)  ae  =  16a*. 


the  Scott-Zhang  quasi-interpolation  operator  in  f \  and  we  let  Qp  :  L2(V)  -A  V~onf(r)  be  the 
corresponding  Scott-Zhang  operator  on  T  C  d flj.  We  now  recall  the  definition  and  the  main 
properties  of  these  interpolators.  For  any  p  6  AT  (7^(fli)),  choose  some1  T  C  ujp.  Let  {Ap^  : 
i  =  1,  •  •  •  ,  d  +  1}  be  the  barycentric  coordinates  of  T  and  let  denote  by  {&T,i  ■  i  =  1,  ■  ■  •  ,  d  +  1} 
its  L2-dual  basis,  i.e. ,  (A =  $i,j-  Let  {<l>p}pejyf T~(fi:))  denote  the  set  of  nodal  basis  of 

V~on{(fii).  Then,  the  Scott-Zhang  quasi-interpolation  operator  is  defined  by 

Qirn  =  (f  dp'n)  ’ 

P&M(T~hm)  v  T 

The  operator  Qp  is  defined  similarly.  Both  operators  enjoy  the  following  approximation  and 
stability  properties: 

Lemma  A. 2.  For  any  g  £  H the  quasi-interpolation  operators  Qi  :  i71(f \)  -A  V~onf(flj) 
and  Qp  :  L2(r)  -A  V~onf(r)  with  T  C  di\,  satisfy  the  following  properties: 

II  Qidl|o,T  <  Iloilo, WT)  IIQp?IIi,t  <  llhlli  ,uT  >  \\Qrv\\o,F  ;$  llhllo  ,oF  •  (A. 2) 

-  QiV)\\o,T  <  ||hl|l,^T  >  II^U  -  Qrv)\\o,F  <  WvWl, Of  •  (A. 3) 

"^Note  that  the  choice  of  T  is  not  unique 
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(a)  \\E\\2As  for  NIPG:  ae  =  4a* 


levels 

h 

e 

10-° 

lO-'4 

1CT3 

urz 

10”1 

1 

101 

1CH 

103 

104 

10b 

0 

2_1 

1.6 

1.6 

1.6 

1.6 

1.6 

1.6 

1.6 

1.6 

1.6 

1.6 

1.6 

1 

2~z 

1.3 

1.3 

1.3 

1.3 

1.3 

1.3 

1.3 

1.3 

1.3 

1.3 

1.3 

2 

2"3 

1.4 

1.4 

1.4 

1.4 

1.4 

1.4 

1.4 

1.4 

1.4 

1.4 

1.4 

3 

2-4 

1.4 

1.4 

1.4 

1.4 

1.4 

1.4 

1.4 

1.4 

1.4 

1.4 

1.4 

(b)  \\E\\as  for  NIPG:  ae  =  8a* 


levels 

h 

e 

10-b 

10“4 

10"3 

io-'z 

10"1 

1 

101 

10z 

103 

104 

10° 

0 

2_1 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1 

2~z 

0.83 

0.83 

0.83 

0.83 

0.83 

0.83 

0.82 

0.83 

0.83 

0.83 

0.83 

2 

2"3 

0.86 

0.86 

0.86 

0.86 

0.86 

0.86 

0.85 

0.86 

0.86 

0.86 

0.86 

3 

2-4 

0.88 

0.88 

0.88 

0.88 

0.88 

0.88 

0.87 

0.88 

0.88 

0.88 

0.88 

(c)  \\E\\2As  for  NIPG:  ae  =  16a* 


levels 

h 

e 

10“° 

10"4 

10"3 

10~'z 

10"1 

1 

101 

10z 

103 

104 

10° 

0 

2-1 

0.68 

0.68 

0.68 

0.68 

0.68 

0.67 

0.68 

0.68 

0.68 

0.68 

0.68 

1 

2~z 

0.60 

0.60 

0.60 

0.60 

0.60 

0.60 

0.60 

0.60 

0.60 

0.60 

0.60 

2 

2"3 

0.61 

0.61 

0.61 

0.61 

0.61 

0.61 

0.61 

0.61 

0.61 

0.61 

0.61 

3 

2-4 

0.61 

0.61 

0.61 

0.61 

0.61 

0.61 

0.61 

0.61 

0.61 

0.61 

0.61 

Table  7.11.  Norm  of  the  error  propagator  E  =  (I  —  Ds  lA)  for  A  corresponding 
to  NIPG  discretization,  with:  (a)  ae  =  4a*,  (b)  ae  =  8a*  and  (c)  ae  =  16a*.  . 


Furthermore,  both  operators  are  linear  preserving;  i.e.  if  77  G  V~oni(ujT),  then  Qip(x)  =  r](x), 
and  similarly  if  77  G  V~oni(Oe),  then  Qrp(x)  =  rj(x). 

We  refer  to  [60,  66]  for  a  proof  of  the  above  result. 

By  using  Q,  and  Qr,  we  can  now  define  our  local  interpolation  operator  P^,  for  which  the 
properties  (5.23)-(5.24)  can  be  shown.  The  definition  is  inspired  on  the  ideas  from  [15]  for  the 
weighted  L2-projection.  To  define  such  an  operator,  we  treat  the  interior  of  each  subdomain  and 
the  interfaces  differently.  We  define  the  interpolation  operator  Pfr  :  V^R  — >  V~oni  as  follows: 

{QiEiV,  at  vertices  inside  the  subdomain  f \ 

QrEiv,  at  vertices  inside  each  face  T  of  the  subdomain  f 2*  (A. 4) 

0,  at  vertices  elsewhere  (cross-points). 

The  next  result  guarantees  that  the  operator  Pfc  defined  above  does  satisfy  the  approximation 
and  stability  properties  (5.23)-(5.24): 

Lemma  A. 3.  For  any  v  G  V^R,  the  operator  Pjf  :  V^R  -A  4~conf  satisfies 

Wi1  ~  pbv\\o,K  <  h\  \og2h/h\l/2\\v\\iXK, 

\Phv\i,K  <  |  log2h//i|1/2||u||i)ft.!K  . 


(A.5) 

(A.6) 
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(a)  Preconditioned  GMRES  for  IIPG:  ae  =  a* 


levels 

h 

10"5 

lO"4 

icr3 

10"2 

10-' 

1 

101 

102 

103 

104 

105 

0 

2-i 

X 

X 

6 

5 

4 

3 

3 

2 

2 

2 

1 

1 

2~2 

X 

6 

5 

4 

3 

2 

2 

2 

2 

2 

2 

2 

2~d 

X 

X 

5 

4 

3 

3 

2 

2 

2 

2 

2 

3 

2-4 

X 

X 

5 

4 

3 

3 

2 

2 

2 

2 

2 

(b)  Preconditioned  GMRES  for  IIPG:  ae  =  2a* 


levels 

h 

10"5 

lO"4 

1CT3 

1CT2 

10-1 

1 

101 

102 

103 

104 

10s 

0 

2_1 

X 

4 

3 

3 

2 

2 

2 

1 

1 

1 

1 

1 

2"2 

X 

4 

3 

2 

2 

2 

2 

1 

1 

1 

1 

2 

2"3 

X 

4 

3 

3 

2 

2 

2 

2 

2 

1 

1 

3 

2-4 

X 

4 

3 

3 

2 

2 

2 

2 

2 

2 

2 

(c)  Preconditioned  GMRES  for  IIPG:  ae  =  4a* 


levels 

h 

10-5 

10-4 

10"3 

10"2 

10-1 

1 

101 

102 

103 

104 

105 

0 

2_i 

4 

3 

3 

2 

2 

1 

1 

1 

1 

1 

1 

1 

2"2 

30 

3 

2 

2 

2 

1 

1 

1 

1 

1 

1 

2 

2"3 

16 

3 

2 

2 

2 

2 

1 

1 

1 

1 

1 

3 

2-4 

11 

3 

2 

2 

2 

2 

1 

1 

1 

1 

1 

Table  7.12.  Number  of  GMRES  iterations  with  the  preconditioner  Asl  for 
IIPG,  with:  (a)  ae  =  a*;  (b)  ae  =  2a*;  and  (c)  ae  =  4a*. 


The  proof  follows  the  ideas  from  [15,  Lemma  4.6]. 


Proof.  We  start  showing  the  approximation  estimate.  Using  standard  triangle  inequality,  the 
stability  and  approximation  properties  of  Ei  given  in  Lemma  A.l,  together  with  the  approxi¬ 
mation  result  (A. 3)  of  the  Qi  from  Lemma  A. 2,  we  have 

\\v  -  p}v Ho^  <  || v  -  QiEiV Ho,^  +  ||  QiEiV  -  pj}v\\ o,^ 

<  ||u  -  EivWo&i  +  || (I  -  Qi)EiV Ho.n.;  +  ||  QiEiV  -  p}v\\ o,^ 

<  h\v\i)h&i  +  h\\Eiv\\itni  +  || QiEiV  -  P^v\\ o,^  , 

<  h\v 1 1,^  +  h\\v\\i,h, a  +  II QiEiV  -  ,  (A. 7) 


Hence,  to  conclude  we  only  need  to  estimate  J|| QiEiV  —  Pftv ||o,ar  To  simplify  the  notation, 
throughout  the  proof  we  set  x  =  G  ^onf  as  defined  in  (A. 4),  and  denote  Xi  :=  QiPiV- 
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(a)  Preconditioned  GMRES  for  NIPG:  ae  =  8a * 


levels 

h 

10"5 

lO"4 

1CT3 

10"2 

10-1 

1 

101 

102 

103 

104 

105 

0 

2_i 

4 

3 

3 

2 

2 

2 

1 

1 

1 

1 

1 

1 

2"2 

X 

3 

3 

2 

2 

2 

1 

1 

1 

1 

1 

2 

2-3 

X 

3 

3 

2 

2 

2 

2 

1 

1 

1 

1 

3 

2-4 

X 

3 

3 

2 

2 

2 

2 

2 

1 

1 

1 

(b)  Preconditioned  GMRES  for  NIPG:  ae  =  16a* 


levels 

h 

10"5 

1CT4 

10"3 

10"2 

10-1 

1 

101 

102 

103 

104 

105 

0 

2_1 

3 

3 

2 

2 

2 

1 

1 

1 

1 

1 

1 

1 

2"2 

5 

3 

2 

2 

2 

1 

1 

1 

1 

1 

1 

2 

2"3 

3 

3 

2 

2 

2 

1 

1 

1 

1 

1 

1 

3 

2-4 

4 

3 

2 

2 

2 

1 

1 

1 

1 

1 

1 

Table  7.13.  Number  of  GMRES  iterations  with  the  preconditioner  Asl  for 
NIPG,  with:  (a)  ae  =  8ct*  and  (b)  ae  =  16a*. 


Then  by  using  discrete  L 2  norm,  we  have 
II  QiEiV  -  P%v ||oa  =  || X  -  XilloA 

<  hd(X-Xi)2{p) 

pedCiinMiiCii 

<  V  Ekd(x-xt)*(p) 

r cdQi per 


%  5Z hd  (Xi  -  QrEiy)2  (p)  +  hdx2i(p) 

redo,  \per  pear 

~  h\\(Xi  ~  QrEjv)\\Qe  +  II Xi llo,ar  j 

redOi  Veer  / 


(A.8) 


(A.9) 


We  need  to  bound  the  two  terms  appearing  in  the  last  expression.  Taking  into  account  that  Qr 
is  linear  preserving  (cf.  Lemma  A. 2)  and  also  using  its  Testability  property  (A. 2),  the  standard 
trace  inequality,  and  the  local  approximation  property  of  the  Scott-Zhang  interpolation  Q,;  (A. 3), 
we  find  for  the  first  term 

h\\Xi  ~  QrEiV\\le  =  h\\QrXi  ~  2r-B^||o,e  Z  ^llx*  “  EM\l,oe 

<  h  (diam(o;p))_1  || Etv  -  QiE^ ||o)Wp  +  hdi&m(u)p)\EiV  -  QiEiV \\^p 
<hdia.m(ojp)\\Eiv\\liUp, 

where  Oe  C  di\  is  the  local  patch  associated  with  e  €  7jl\dni  on  the  interface,  and  up  C  fij  is 

the  patch  associated  with  the  vertex  p£e.  Note  that  diairi(wp)  =  2 h. 

Summing  up  the  above  inequality  for  all  edges/faces,  we  obtain  that 

^  X^llx*  -  Qr^«||g,e  Z  Miam(£Up)||E^||?A  <  hdiam(ujp)\\v\\lAn..  (A.10) 

rcdRi  ecr 


38 


B.  AYUSO  DE  DIOS,  M.  HOLST,  Y.  ZHU,  AND  L.  ZIKATANOV 


To  bound  the  second  term  in  (A. 9)  we  have  to  distinguish  between  the  2D  and  3D  cases.  In  the 
2D  case,  T  is  one-dimensional,  (it  is  a  set  of  “edges”)  and  so  9T  reduces  to  its  two  endpoints, 
say  {p,  q}.  Hence, 

IIXi|lo,ar  =  (IXi(p)l2  +  \Xi(q)\2)  <  llx*llo,oo,u!p  +  ||Xi|lo,oo,«,  >  5r  =  {P>9>  » 

To  bound  each  of  the  above  two  terms  on  the  right  side,  we  use  the  two-dimensional  discrete 
Sobolev  inequality  [15,  Lemma  2.3]; 


||Xi||o,oo,ajp  E  c  log 


diam(o;p)  \ 


1 1 X*  1 1  l,c* 


(A. 11) 


and  so  summing  over  all  T  G  the  resulting  estimate,  and  using  the  stability  (A. 2)  from 
Lemma  A. 2  of  Q%  together  with  the  stability  and  approximation  properties  of  E{  given  in 
Lemma  A.l,  we  finally  get 


E  HXillo,ar  £  E  E  log 

r&dUi  reafi; pear 


llxillu^  <log  —  ||x»lli,n,- 


diam(cup 


=  log  —  \\QiEiv\\ltQi  <  log  —  [|u||1A  . 


(A.  12) 


Next,  we  give  the  corresponding  estimate  for  the  3-dimensional  case.  In  3D,  any  T  C  is  a 
set  of  faces  and  edges.  Notice  however,  that  it  is  enough  to  consider  the  case  when  T  is  two- 
dimensional  (face),  since  the  other  case  reduces  to  the  estimate  already  done  for  2D.  T  being 
a  “face”,  <9T  reduces  to  a  set  of  edges,  {e  :  e  C  <9T}.  Hence,  defining  the  set  coe  :=  {T  G  7~  : 
Tne/0eC<9T}  and  using  now  the  discrete  Sobolev  inequality  [15,  Lemma  2.4]  (instead  of 
(A. 11)),  we  find 

llx.llo,ar=  E  ItallLS  E  log  (dla^«))  ||Xi||jw  . 

ecar  ecar  ^  ' 

Summing  the  above  estimate  over  all  T  G  <912*  and  using,  as  before,  the  stability  (A. 2)  from 
Lemma  A. 2  of  Q-L  together  with  the  stability  and  approximation  properties  of  E)  given  in  Lemma 
A.l,  we  find, 


X]  llx<llo,ar  £  E  E  log 


redtti ecar 


diam(we) 


IIXilll,„.  S  log 


E  E  ii villa 


reaiti  ecar 


;$ log  y~h  J  ~ log  y~h  J  • 

Now,  substituting  into  (A. 9)  the  above  estimate  (or  correspondingly  (A.  12))  together  with 
(A.  10)  we  finally  get 


II QiEiV  -  Phv ||0)Qi  —  || X  -  Xillo.n,  —  C2Ch  |H|1;/iq.  +  Ch  log  J  |M|1;f2.  . 

Inserting  this  estimate  in  (A. 7),  the  proof  of  the  approximation  property  is  concluded  by  using 
now  the  definition  of  the  il1-weighted  norm. 

Finally  we  show  the  stability  of  P^  (A. 6).  From  the  definition  of  the  norm 


I  P^SjII2  ^ 
I  -^h  ^  II  1,k,£1 


E 


hv  111, Hi  1  \\Phv\\l,a 


E 

Tcf2,  ,|J  T=Qi 


ph„ ,  1 1 2 
^hv\\l,T 
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Note  that  Pftv  G  and  v  G  V^R.  To  deal  with  possibly  different  mesh  sizes  we  consider 

the  local  L2-projection  Vt  '■  L2(T)  — >  P1(T)  for  any  T  G  7r.  For  h  >  h,  such  an  element  is 
the  union  of  other  subelements  in  the  partition  Th-  Then,  adding  and  subtracting  Vtv,  triangle 
inequality  together  with  inverse  inequality  and  the  approximation  property  (A. 5),  gives 

\Phv\i,T  <  | Pjfv  —  Vtv\i,t  +  \Ptv\i,t  <  C'(Zi)-1 \\Pftv  —  Vtv\\o,t  +  \Ptv\i,t 


<  C(h )  (J| Phv  -  u||o,t  +  lb  -  Ptv\\o,tj  +  C\v\1:T  <  C(h)  || Pftv  —  i>||o,r  +  C'IM|i,t  • 

The  Stability  now  follows  immediately,  by  summing  over  all  elements  T  C  b,  using  the  definition 
of  the  weighted  H 1  -semi-norm  and  the  weighted  L2-norm  together  with  the  approximation  result 
already  shown: 

\Phv\i,K,n  <  Ch~l\\p};v  -  u||0,«,n  +  lb||i,feA,n  <  Ch~lh  f^log  IMIiAaP  +  IMIiA^n 


b  log 


and  the  proof  is  complete. 


Appendix  B.  Proof  of  (5.19)  and  (5.18) 

Since  the  results  in  this  section  concern  only  the  space  VfR,  we  will  omit  the  superscript  CR 
and  write  Ao  instead  of  Aq  ,  <pe  instead  of  (p^R ,  etc.  We  also  point  out  that  the  two  Lemmas 
that  follow  are  cases  of  much  more  general  theorem  on  additive  methods  which  can  be  found  in 
many  texts  (see  e.g.  [77],  [69],  [67]). 

We  first  prove  the  identity  for  the  Jacobi  method  (5.19). 

Lemma  B.l.  Let  R~x  be  defined  via  the  expression 


R~1Aqv  =  J2 


Ao(v,ipe) 
A0(<p  e?  <Pe) 


Then  for  all  w  G  V)yR  the  following  relation  holds: 

(Rw,w)  =  ^  Ao(<Pe,<Pe 


Proof.  Given  w  G  VffR,  by  setting  X  =  [R  1Aq\  1  we  obtain  that 


w  =  R  1Aq[R  1  Aq]  lw  =  R  1Aq(Xw)  =  ^22 


'  Ao(X'W,(pe) 
O  A0{(Pe,Te) 


In  other  words,  we  have  w  =  Y!e&£°weTei  with  we  =  2222’’^)  '  then  easily  obtain  the 
desired  result  by  the  following,  rather  straightforward,  computations: 

A  /  ,  X  .2  [Ao{Xw,  ipe)}2  Ao{Xw,ipe) 


^  A0(pe,Pe)l 


1o  Ao{<Pe,ipe) 


y  weAo(Xw,<pe)  =  Ao(Xw,  y  we<pe)  =  Ao(Xw,w) 


=  (Aq[R  1Aq]  Lw,w)  =  ( Rw,w ). 
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The  following  Lemma  verifies  the  identity  for  ( B  lv,v)  given  in  (5.18). 
Lemma  B.2.  Let  B  be  defined  as  in  (5.13).  We  then  have 

( B~lv ,  v )  =  inf  [K(v  ~x,v~x)  +  a(x,  x)]- 

XSY-Conf 


Proof.  Let  x  £  L~onf  be  arbitrary  and  rj  =  {Ac)  lQc B  xv.  It  is  clear  that  Acr]  =  Qc B  lv 
and  moreover, 

v-ri  =  v-  ( Ac)~1QcB~1v  =  BB-'v  -  PcA-lB~lv  =  (B  —  ( Ac)~1Qc)B~1v  =  R-lB~lv. 
Setting  i] o  =  X  ~  7h  so  that  X  =  V  +  ??o  then  shows  that 

K(v-x,v~x)  =  (R(v-r]-r]o),v-r)-r]o) 

=  (R-^RB-'v-Tio^iRB-'v-Tio)) 

=  (RB~lv,  B~lv)  -  2(770,  B~lv)  +  {R~lri 0,  %), 

and  also 

a(x,x)  =  a(v  +  Vo),V  +  V  0) 

=  a(v,  V)  +  2 a(7/,  ??0)  +  a(??0,  %) 

=  (^G'7,  ??)  +  2(Acii,  r] 0)  +  0(7/0, 770) 

=  (QCB~1v,  (Ac)~1Qcri)  +  2(B-lv,ixf)  +  0(7/0, 7/0). 

Adding  the  last  two  identities  and  using  the  fact  that  0(7/0,770)  >  0  and  (i??/ 0,770)  >  0,  and 
applying  the  definitions  of  anc[  77  (in  that  order)  then  gives 

7l(v  -  x,v  -  x)  +  a(x,x)  >  {RB~lv,B~1v)  +  {QcB~1v,{Ac)~lQcB~1v) 

=  {B-'v,  RB~lv)  +  {B~lv,  ( Ac)~lQcB~1v ) 

=  {B-'v,  (R  +  ( Ac)~1Qc)B~1v )  =  (B_1u,  u) 

v - V - ' 

B 

The  proof  is  complete  because  x  was  arbitrary  and  moreover,  equality  holds  if  and  only  if 
7/0  =  0.  □ 
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